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Foreword 


This program is an important tool in the study of alternative routes between the Earth and 
the Moon. Dr. John Aired was the NASA Technical Monitor for the contract under which this 
program was produced. Mr. Andy Petro was the NASA Task Manager for this particular task. 
Mr. Bill Stump was the Eagle Project Manager for the contract under which this program was 
produced. The program was written by Jack Funk, originally in Quick Basic, and translated into 
Fortran by Mr. Bill Engblom. Mr. Engblom also prepared the documentation. 
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1.0 Introduction 


The program LP1 calculates outbound and return trajectories between low earth orbit 
(LEO) and libration point #1 (LI). Libration points (LP) are defined as locations in space that 
orbit the Earth such that they are always stationary with respect to the Earth-Moon line. LI is 
located behind the Moon such that the pull of the Earth and Moon together just cancel the 
centrifugal acceleration associated with the libration point’s orbit. 

The outbound flights depart from a circular orbit of any altitude and inclination about the 
Earth and culminate in a circular orbit about the Earth at libration point #1 within a specified 
flight time. The flight involves three bums. 

First, the departure orbit is made into a more eccentric orbit (ellipse or hyperbola) with 
an initial AV in order to reach the lunar sphere of influence (SOI), a region where the vehicle is 
near its lowest velocity in the trajectory. The SOI is a spherical region whose surface normally 
includes all the points at a distance of 11% of the Earth-Moon distance from the Moon’s center. 
However, in order to simplify the calculations this radius was increased to include LI, enlarging 
the sphere radius to 15% of the Earth-Moon distance. Note, this change in SOI radius should not 
change the results significantly. A given flight may penetrate the SOI at a number of points 
identified by projecting lunar latitude and longitude onto the SOI. For each flight the program 
will calculate a set of possible trajectories associated with a set of SOI penetration points — a 
matrix of longitudes and latitudes. 
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Next, a second bum (at the SOI) is performed involving a "flyby" of the Moon from the 
SOI point above the front side of the Moon to LI behind the back side. Note, the SOI penetra¬ 
tion point and LI will always be the same distance from the center of the Moon. From the 
geometry of the trajectory it is apparent that the Lunar "flyby" perigee altitude, supplied by the 
user, will occur midway between the SOI point and LI. Consequently, the geometry of the orbit 
will force true anomaly, flight path angle, and absolute time to or from perigee passage to be the 
same for both the SOI and LI points. There are two paths between these points, posigrade and 
retrograde. LP1 calculates only posigrade "flyby" trajectories since retrograde orbits that pass 
through the perigee altitude are not always possible while there is always a posigrade solution. 
The retrograde trajectories warrant more study in future work. Since LI is constantly rotating 
with the Moon, this trajectory is iterated until LI is reached. 

The third bum is simply a circularization of the trajectory at LI about the Earth. 

The velocity vector is corrected to that of LI. Note, once the SOI to LI trajectory has been 
established, the Earth-SOI flight is iterated until the total transfer time, including the transfers 
from LEO to the SOI, and SOI to LI, match the user’s flight time constraint (an input value). 
This is done for a matrix of SOI penetration points, as mentioned earlier. 

The return trajectories, which start at LI and finish in the specified LEO orbit within the 
specified flight time, are calculated similarly. For instance, the "flyby" trajectory is calculated 
first, starting at LI and finishing at the SOI via a posigrade orbit calculated using the same 
geometric simplifications described above. Note, the flight profile used in this program to 
calculate flights between Earth and LI may not be the optimum with respect to a minimum AV. 
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After the user has defined the trajectory as outbound or return, the Earth orbit altitude 
and inclination, and the total flight time, LP1 produces matrices which display the total AVs, the 
three component AVs (described above), the "fly-by" trajectory inclinations, and the "fly-by" 
azimuth angle at the SOI for the resulting flight from Earth to LI for a representative set of SOI 
points. These points are defined by the user, who provides a starting longitude and latitude, and 
an increment for each. The matrix is built with 10 longitudes forming the columns and 19 
latitudes forming the rows. 

Section 2.0 of this document describes the input required from the user to define the 
flight. Section 3.0 describes the contents of the six reports that are produced as outputs. Section 
4.0 includes the instructions needed to execute the program. 

A more detailed description of the process used in LP1 has been included as Appendix D 
(main program). Appendix E (in-program subroutine), and Appendices F, G, H, I, and J (external 
subroutines). LP1 was derived from the PLANECHG program (also produced under this 
contract) with the major addition of the FLYBY subroutine. Therefore, the documentation for 
PLANECHG may be used as a reference for many of the equations, variables, and conventions 
used in LP1 (except in the FLYBY routine). 
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2.0 Program Inputs 


The following paragraphs discuss the inputs provided by the user. The prompt is the 
message displayed by the program onto the screen. The input variable is the program variable 
assigned to the user’s response. The description provides information about how to respond to 
the prompt. 

1. Prompt: INPUT OUTBOUND OR RETURN 

Input variable: MD 

Description: Enter OUTBOUND for Earth-to-Ll trajectory calculations. Enter 
RETURN for Ll-to-Earth trajectory calculations. 

2. Prompt: INPUT PERIGEE ALTITUDE OF EARTH ORBIT (NMI) 

Input variable: HPE 

Description: Enter the height above the Earth’s surface of the Earth circular orbit, in 
nautical miles. 

3. Prompt: INPUT PERIGEE ALTITUDE OF LUNAR ORBIT (NMI) 

Input variable: HPM 

Description: Enter the height above the Lunar surface of the Lunar circular orbit, in 


nautical miles. 




4. Prompt: INPUT EARTH DEPARTURE JULIAN DATE 

Input variable: TIMJ 

Description: Enter the origin trans-SOI injection date (Earth departure date for 
outbound trajectories) in Julian day format, where January 1, 2000 is day 
2,451,545. Refer to Section C of "The Astronomical Almanac of the Year 
1988". Day 2451545 is the default value if zero is entered in this field. 

5. Prompt: INITIAL LONGITUDE 

Input variable: ALONI 

Description: Enter initial sphere of influence longitude for the output matrices. This 
value will become the heading for column 1 of the matrices. 

6. Prompt: INPUT INCREMENT FOR THE MAP 

Input variable: DELLON 

Description: Enter longitude increments for the output matrices. Applied to the initial 
longitude, this value defines the subsequent column headings of the 
matrices. Longitudes for outbound trajectories should be between zero 
and -90 degrees. 
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7. Prompt: 


INPUT INITIAL LATITUDE 


Input variable: ALATI 

Description: Enter initial sphere of influence latitude for the output matrices. This 
value will become the heading for row 1 of the matrices. 

8. Prompt: INPUT INCREMENT FOR THE MAP 

Input variable: DELLAT 

Description: Enter latitude increments for the output matrices. Applied to the initial 
latitude, this value defines the subsequent row headings of the matrices. 


9. Prompt: INPUT EARTH ORBIT TO LUNAR ORBIT INCLINATION 

Input variable: AINCEO 

Description: Enter Earth circular orbit inclination, in degrees. This is not what is 
typically considered inclination (i.e., a measurement taken from the 
Earth’s equatorial plane), but rather the angle between the plane of the low 
Earth orbit and the plane of the Moon’s orbit about the Earth. 

10. Prompt: INPUT FLIGHT TIME 

Input variable: FTIM 

Description: Enter the desired total flight time from LEO to LI, in hours. 
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3.0 Program Outputs 


This section describes the contents of each of the six reports generated by the program. 
These reports may be found in the output file, LPl.OUT. Samples of all six reports have also 
been included. 

Report #1 : Total Delta Velocity Map For Outbound/Retum Trajectories 

The top section of this report repeats the input values entered by the user. The second 
section is a 10 X 19 matrix of total velocities (in ft/sec) required to fly the profile described by 
the inputs. Each cell corresponds to a particular latitude and longitude on the Lunar sphere of 
influence. A "flyby" bum (transfer from SOI to LI) may occur at any one of these coordinates, 
and the value of the corresponding cell is the total velocity required for the flight if the "flyby" 
bum occurs at that location. The total is the sum of the Earth-SOI transfer orbit injection AV 
(from LEO to SOI), the sphere of influence to LI, Lunar "flyby" AV (from SOI to LI), and the 
destination circular orbit injection AV (circularization at LI). Note, the longitudes must always 
be between 0° and -90° for outbound trajectories, and between 0° and +90° for return trajec¬ 
tories. The third section of the report is a summary of key data corresponding to the matrix cell 
containing the lowest total velocity. This data includes: 

• X-, Y-, and Z-components of velocity at the sphere of influence just before and just 

after the "flyby" bum (VX, VY, VZ). 

• Total magnitude of the velocity at the sphere of influence just before and just after the 
"flyby" bum (VEL). 

• Flight path angle at the sphere of influence just before and just after the "flyby" bum 

(GAMA). 

Azimuth of the sphere of influence point from the Earth and from the Moon (AZM). 
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• Earth and "flyby" orbit inclinations measured from the Earth-Moon plane (AINC). 

• Earth-SOI transfer trajectory and "flyby" orbit insertion ascending/descending node 
positions measured from Earth-Moon line at 0° longitude (ANODE). 

• Earth-to-SOI and SOI-to-Ll times of flight (TIME). 

■ Earth-SOI transfer trajectory and final LI orbit insertion AV’s (DVCIR). 

• Sphere of influence, Lunar "flyby" AV (DVPHER). 

• Total AV (DVTOTAL). 
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Re port #2 

Map of Delta Velocity at Sphere of Influence for Outbound/Retum Trajectories 

This report is a matrix of the delta velocities that occur at the sphere of influence to 
match the velocity vectors of the Earth-to SOI trajectory and the SOI-to-Ll "flyby" trajectory for 
each SOI point. Each cell corresponds to a particular latitude and longitude on the sphere of 
influence at which a bum may occur. 
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DELTA VELOCITY AT SPHERE OF INFLUENCE FOR EARTH TO LI FLYBY TRANSFER TRAJECTORIES NODE 
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Report #3 

Map of Inclinations of Flyby Orbits for Outbound/Retum Trajectories 

This report contains a matrix of the inclinations of the "flyby" trajectories with respect to 
the Earth-Moon plane, between each SOI point and LI. Each cell corresponds to a particular 
latitude and longitude on the sphere of influence at which the trajectory originates or culminates. 
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OF INCLINATION OF FLY-BY ORBIT FOR EARTH TO LI FLYBY TRANSFER TRAJECTORIES NODE OPTION 
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Report #4 

Map of Delta Velocity at LI for Outbound/Retum Trajectories 

This report is a matrix of the delta velocities that occur at LI to circularize the trajectory 
at LI (outbound flights), or to initiate a posigrade "flyby" trajectory from LI to a particular SOI 
point (return flights). Each cell corresponds to a particular latitude and longitude on the sphere 
of influence through which the "flyby" trajectory originates or culminates. 
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OF DELTA VELOCITY AT LI FOR EARTH TO LI FLYBY TRANSFER TRAJECTORIES NODE OPTION 
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Re port #5 

Map of Delta Velocity at Earth for Outbound/Retum Trajectories 

This report is a matrix of the delta velocities that occur to initiate a transfer orbit from 
LEO to the SOI. Each cell corresponds to a particular latitude and longitude on the sphere of 
influence at which a bum may occur. 
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DELTA VELOCITY AT EARTH FOR EARTH TO LI FLYBY TRANSFER TRAJECTORY NODE OPTION 0 
DATE 7-NOV-88 TIME 14:46:12 

EARTH TO LI FLYBY TRANSFER FLIGHT TIME (HR) = 200.0 INCL EARTH 25.0 INCL MOON = 90. 


O'fWinNHNUlOOffiTOIOMHNIflCn^O 

rHOOOOOOOOOOOOOOOOOrH 

ooooooooooooooooooo 


OC'J*a , t''T-IV£>CM<T>COCOCOO'><N'£>pir-' 3 , CMO 

oooor~r-ou>ininininin'o'nr-t-~oooo 

rHOOOOOOOOOOOOOOOOOrH 

ooooooooooooooooooo 


OOOOCM^OOni-tOrHCOOO^rO^OOOO 

rHOOOOOOOOOOOOOOOOOrH 

ooooooooooooooooooo 


o'CfnoffioooooTn^oonomon^o 

ocor~vo"5r'<j>mcMCMCMcs!CMcn''* , 3'vor-ooo 

rHOOOOOOOOOOOOOOOOOrH 

OOOOOOOOOOOOOOOOOOO 

HrHHHHrHrHHHHHrHrlHrHHrlHH 

O’Tor'tO'-ooooooocnvoinr^O'q'o 
oaor'in^rrocM cm m rr in r~ ao o 

rHOOOOOO OOOOOOrH 

ooooooo ooooooo 

rHrHi-HrHrlrHrH rHrHHrlrlrHrH 

oror-'srro^ooooooo^ro^pr'mo 
o oo io m n nv in ® oo o 

rHOOOOO OOOOOrH 

OOOOOO OOOOOO 


OCMinCHrHOOOOOOOOOrHCMinCMO 

o oo o in in io oo o 

rHOOOO OOOOrH 

ooooo ooooo 

i—I r— I rH i —I t— } 1-HrHrHrHrH 

OrH^OOOOOOOOOOOOO^rHO 

oooioin^r rr m ® oo o 

rHOOOO OOOOrH 

ooooo ooooo 


ooooooooooooooooocooo 
m oo id o oo o 

O O O O O rH 

o o o o o o 


a 


ooooooooooooooooooo 


ooooooooooooooooooo 
ooor-oin'a'mcMrH rHcsiro^inmr-ooo 

I I I I I I I I I 














Report #6 

Map of Azimuths of Flyby Orbits at SOI for Outbound/Retum Trajectories 

This report contains a matrix of the azimuth angles of the "flyby" trajectories at the SOI 
for each particular SOI point. Each cell corresponds to a particular latitude and longitude on the 
sphere of influence through which the "flyby" trajectory originates or culminates. 




OF AZMM FOR EARTH TO LI FLYBY TRANSFERTRAJECTORY NODE OPTION 0 
DATE 7-NOV-8 8 TIME 14:46:12 
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4.0 Program Execution Instructions 


The following instructions describe the steps to be taken by the user to execute this program: 

A. Obtain access to the DEC VAX minicomputer and sign on with user identification. 

B. At the $ prompt, type RUN LP1 . 

C. When prompted by the program, enter the program inputs. See section 2.0 for a 
discussion of the inputs. 

D. After the last input has been entered, the program will execute for approximately 5 
minutes, during which comments will appear notifying the user which trajectory, by 
SOI latitude and longitude, is currently being calculated. Upon completion, the 
'message FORTRAN STOP will appear, followed by the $ prompt. 

E. The program outputs will be placed in a file named LIB RATE. O UT ;### where ### is 
a system generated number of the report. To print the most recently generated report, 
type the following at the $ prompt: TYPE LPl.OUT 

F. To re-execute the program with new parameters, begin again at step (B) above. 
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Appendix A. Program Flow Chart 
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Appendix B. Code Listing 




c ************************************************************ 
C *** Libration Point Program (between LI and Earth) ** 

Q ************************************************************ 

C * Written in Quick Basic by: Jack Funk 
C * Translated by Bill Engblom 

C * Documented by Bill Engblom 

q ************************************************************ 

C 

C MAIN 

C 

IMPLICIT REAL*16 (A-H,0-Z) 

CHARACTER*10 MD, TRAJ, TRAJM, TRAJE 
CHARACTER*8 TIMP 
CHARACTER*9 DATP 
CHARACTER*26 HEAD 

DIMENSION DELV(10,19), VELMOUT(10,19),ALONO(10),ALATP(19) 

* PAGEl(lO), PAGE2(10,19 ), PAGE3(10,19),PAGE4(10,19), 

* PAGE5(10,19), PAGE6(10,19) 

C 

C OPEN OUTPUT FILE 

C 

OPEN (UNIT = 1 f FILE = 'LPl.OUT' , STATUS = '’NEW'’) 

C OUTPUT TO FILE; IP: OUTPUT TO SCREEN; IS 

IP ■ 1 
IS =5 

DPR - 57.29578 
PI = 3.1415926535 
CMUE = 1.407647E+16 
CMUM = 1.731432E+14 
FTNM » 6076.115 
REE = 20925741. 

REMO = 5.7039E+6 
RREM = 207559. * FTNM 
FTM = .3048 

PRINT *, 'INPUT OUTBOUND OR RETURN ' 

READ (6, 5) MD 
5 FORMAT(A10) 

IF (MD .EQ. 'RETURN' )THEN 

MOOD = -1 

HEAD = 'Ll TO EARTH FLYBY TRANSFER' 

ELSE 

MD = 'OUTBOUND' 

MOOD = 1 

HEAD = 'EARTH TO LI FLYBY TRANSFER' 

END IF 

C PRINT *,'INPUT NODE OPTION 1 OR 2 ' 

C READ *, NP 

10 CONTINUE 

PRINT *, 'INPUT PERIGEE ALTITUDE OF EARTH ORBIT (NMI) ' 






READ *, HPE 

PRINT *, 'INPUT PERIGEE ALTITUDE OF LUNAR ORBIT (NMI) ' 
READ *, HPM 

RPE = HPE * FTNM + REE 

RPM = HPM * FTNM + REMO 

PRINT *, 'INPUT EARTH DEPARTURE JULIAN DATE ' 

READ *, TIMJ 

IF (TIMJ .EQ. 0.) TIMJ = 2451545. 

PRINT *, 'LONGITUDE FOR OUTBOUND TRAJECTORIES SHOULD BE 

* BETWEEN 0 AND -90 DEG' 

PRINT *, 'AND RETURN TRAJECTORIES BETWEEN 0 AND +90 DEG' 
PRINT *, 'INITIAL LONGITUDE' 

READ *, ALONI 

PRINT *, 'INPUT INCREMENT FOR MAP ' 

READ *, DELLON 

PRINT *, 'INPUT INITIAL LATITUDE' 

READ *, ALATI 

PRINT *, 'INPUT INCREMENT FOR MAP ' 

READ *, DELLAT 

PRINT *, 'INPUT EARTH ORBIT TO LUNAR ORBIT INCLINATION ' 
READ *, AINCEO 
AINCEO = AINCEO / DPR 
AINCE = AINCEO 

PRINT *, 'INPUT FLIGHT TIME ' 

READ *, FTIM 
CALL DATE(DATP) 

CALL TIME(TIMP) 

15 CONTINUE 

WRITE (IP,7) HEAD, NP 

7 FORMAT (T12,' VELOCITY MAP FOR ',A26,' TRAJECTORIES 

* NODE OPTION ',11) 



WRITE 

(IP,17) DATP, TIMP 


17 

FORMAT 

(T27,'DATE',A10, ' TIME',A10) 



WRITE 

(IP,27) TIMJ, HPE, HPM 


27 

FORMAT 

(/,' JUILIAN DAY ', F8.0,' PERIGEE ALT 

(NMI) EARTH 


*F4.0, 

' MOON = ', F6.0) 



WRITE 

(IP,37) FTIM, AINCEO * DPR 


37 

FORMAT 

(' TRANSLUNAR FLIGHT TIME (HR) = ',F5.1, 

' INCL 


* EARTH 

' F5.1,' INCL MOON PAGE3 ') 



IPRINT 

= 0 



II 

= 1 



NN 

= 1 



DVMIN 

= 99999. 



YDLO 

= QSQRT (CMUE / RREM) 



ALAT 

= ALATI / DPR 



ALON 

= ALONI / DPR 



VELM 

= VELMI 



DVSIM 

= 99999. 



21 CONTINUE 





C VEL = VELM 

DELV(NN,II) = 99999. 

PAGEl(NN) =0.0 
PAGE2(NN,II) =0.0 
PAGE3(NN,II) =0.0 
PAGE4(NN,II) =0.0 
PAGE5(NN,II) =0.0 
PAGE6(NN,II) =0.0 
C PRINT *,' ' 

C PRINT ALAT ALON VEL DVT DVSI MOON EARTH 

C *TIMEE TIMEM RREM' 

22 CONTINUE 
C 

C CALCULATE FLYBY TRAJECTORY CHARACTERISTICS 

CALL FLYBY (ALON, ALAT, RREM, RPM, GAMAM, VELM, VPM, AINCM, DELVLP1, 
*AZMM, AVM, TIMM, MOOD) 

ICALL = 1 
GOTO 25 

1000 CONTINUE 

C WRITE (IS,47) ALAT*DPR,ALON*DPR,VELM,DVTOTAL,DELVEL, 

C * TRAJM$, TRAJE$, TIEM, TIMM, RREMER, PHIE*DPR 

C 47 FORMAT (IX,F5.0,2X,F5.0,2X,F6.0,2X,F7.0,2X,F6.0,2X, 

C * A10,2X,A10,2X,F6.1,2X,F6.1,2X,F7.3,2X,F5.0) 

IF (QABS(TIEM) .LE. .0000000000000001) GOTO 23 
C STORE VALUES FOR OUTPUT 

IF(DVTOTAL .LT. DELV(NN,II) ) THEN 
DELV (NN,II) = DVTOTAL 

PAGE1 (NN) = DVTOTAL 

PAGE2 (NN,II) = DELVEL 

VELMOUT (NN,II) = VELM 
ALATP (II) = ALAT * DPR 

PAGE3 (NN,II) = AINCM * DPR 

PAGE4 (NN,II) = DELVLP1 

PAGE5 (NN,II) = DVCIRE 

PAGE6 (NN,II) = AZMM * DPR 

C WRITE (IS,57) ALAT * DPR, ALON * DPR, VELMOUT (NN,II) 

C * ,DELV(NN,II), PAGE2(NN,II), TRAJM, TRAJE, TIEM, 

C * TIMM, RREMER 

C 57 FORMAT (IX,F6.1,IX,F6.1,2X,F7.0,IX,F8.1,IX,F7.1,IX, 

C * A5,IX,A5,1X,F6.1,IX,F6.1,1X,F8.3) 

C WRITE (IS, 67) VXE, VXM, VYE, VYM, VZE, VZM 

C 67 FORMAT (' VXE ',F7.1,' VXM ',F7.1,' VYE ',F7.1,' VYM ', 

C * F7.1, ' VZE ', F7.1, ' VZM ',F7.1) 

END IF 

23 CONTINUE 

ALONO(NN) = ALON * DPR 

IF (NN .EQ. 10 .AND. IPRINT .EQ. 0 ) THEN 

WRITE (IP,77) ALONO(l), ALONO(2), ALONO(3), ALONO(4), 




* AL0N0(5), ALONO(6), ALONO(7),ALONO(8),ALONO(9),ALONO(10) 

77 FORMAT (/,' ALON> ' , 10(2X,F5.1 )) 

WRITE (IP,87) 

87 FORMAT (' ALAT') 

END IF 

IF (NN .LT. 10) THEN 

ALON = (ALONI + DELLON * QFLOAT (NN ) ) / DPR 
NN = NN + 1 
GOTO 21 
END IF 
IPRINT = 1 

WRITE (IP,97) ALAT * DPR, PAGE1(1),PAGE1(2), 

*PAGE1(3),PAGE1(4),PAGE1(5),PAGE1(6), 

*PAGE1(7),PAGE1(8),PAGE1(9),PAGE1(10) 

97 FORMAT (IX, F5.1, 2X, 10(1X,F6.0)) 

IF (II .LT. 19) THEN 

ALAT = (ALATI + DELLAT * QFLOAT (II ) ) / DPR 

II = II+l 
NN =1 

ALON = ALONI / DPR 
GOTO 21 
END IF 

VELM = VELMMIN 
ALAT = ALATMIN 
ALON = ALONMIN 
ICALL = 2 
GOTO 25 
2000 CONTINUE 

WRITE (IP,107) NP 

107 FORMAT (///,T23,'MINIMUM VELOCITY PRINT NODE OPTION ',11) 

WRITE (IP,117) DATP, TIMP 
117 FORMAT (T2 7,'DATE ' ,Al0, ' TIME',A10) 

WRITE (IP, 127) FTIM, ALATMIN * DPR, ALONMIN * DPR, AINCMP * DPR, 

* AINCEO * DPR 

127 FORMAT (/,' FLIGHT TIME = ',F4.0,' LAT = ',F6.1,' LON = ',F6.1, 

* ' LUNAR AINC = ', F5.1,' EARTH AINC = ',F5.1) 

WRITE (IP,137) 

137 FORMAT (' BODY VX VY VZ VEL GAMA AZM 

*AINC ANODE TIME DVCIR' ) 

WRITE (IP,147) VXMP,VYMP,VZMP,VMAGMP,GAMAMP * DPR, 

*AZMMP * DPR, AINCMP * DPR, ANODEMP*DPR,TIMMP,DELVLP1P,DELVELP 
147 FORMAT (' MOON ',4(IX,F6.0),IX,F5.1,IX,F6.1,2X,F5.1,2X,F6.1, 

* IX,F5.1,IX,F7. 0, ' DVPHER — ' ,F7.0) 

WRITE (IP,1475) VXEP, VYEP, VZEP, VELEP,GAMAEP*DPR,AZMEP*DPR, 

* AINCEP * DPR,ANODEEP*DPR,TIEMP,DVCIREP,DVTOTALP 
1475 FORMAT (' EARTH',4(IX,F6.0),IX,F5.1,IX,F6.1,2X,F5.1,2X,F6.1, 

* IX,F5.1,IX,F7. 0, ' DVTOTAL = ', F7.0) 

WRITE (IP,157) CHAR(12) 

157 FORMAT (' ' ,A1) 




WRITE (IP, 167) HEAD, NP 

167 FORMAT (T10,'MAP OF DELTA VELOCITY AT SPHERE OF INFLUENCE 

*FOR ',A26,' TRAJECTORIES NODE OPTION ',11) 

WRITE (IP, 177) DATP, TIMP 
177 FORMAT (T27,'DATE', A10,' TIME',A10) 

WRITE (IP,187) HEAD,FTIM, AINCEO * DPR, AINCM * DPR 
187 FORMAT (' ',A26,' FLIGHT TIME (HR) = ',F6.1,' INCL EARTH ', 

*F6.1,' INCL MOON = ',F5.1) 

WRITE (IP,197) ALONO(l), ALONO(2), ALONO(3), ALONO(4), 

*ALONO(5),ALONO(6), ALONO(7), ALONO(8), ALONO(9), ALONO(IO) 
197 FORMAT (/,' ALON> ',10(2X, F5.1)) 

18 CONTINUE 

WRITE (IP,207) 

207 FORMAT (' ALAT') 

DO 28 NPI = 1 , 19 

WRITE (IP,217) ALATP(NPI), PAGE2(1,NPI),PAGE2(2,NPI), 

* PAGE2(3,NPI),PAGE2(4,NPI),PAGE2(5,NPI),PAGE2(6,NPI), 

* PAGE2(7,NPI),PAGE2(8,NPI),PAGE2(9,NPI),PAGE2(10,NPI) 

217 FORMAT(IX,F5.1,IX,10(IX,F6.1)) 

28 CONTINUE 

WRITE (IP, 227 ) CHAR (12) 

227 FORMAT (' ', Al) 

WRITE (IP, 237 ) HEAD , NP 

237 FORMAT (T6,'MAP OF INCLINATION OF FLY-BY ORBIT FOR ',A26, 

* ' TRAJECTORIES NODE OPTION ',11) 

WRITE (IP, 247 ) DATP, TIMP 

247 FORMAT (T27,'DATE',A10,' TIME',A10) 

WRITE (IP, 257 ) HEAD, FTIM, AINCEO * DPR, AINCM * DPR 
257 FORMAT (' ',A26,' FLIGHT TIME (HR) = ',F6.1,' INCL EARTH 

*F6.1,' INCL MOON = ',F6.1) 

WRITE (IP,267) ALONO(1), ALONO(2), ALONO(3), ALONO(4), 

*ALONO(5),ALONO(6), ALONO(7), ALONO(8), ALONO(9), ALONO(10) 
267 FORMAT (/,' ALON> ', 10(2X, F5.1)) 

WRITE (IP,277) 

277 FORMAT (' ALAT') 

DO 48 NPI = 1 , 19 

WRITE (IP,287) ALATP(NPI),PAGE3(1,NPI),PAGE3(2,NPI), 

* PAGE3(3,NPI),PAGE3(4,NPI),PAGE3(5,NPI),PAGE3(6,NPI), 

* PAGE3(7,NPI), PAGE3(8,NPI),PAGE3(9,NPI),PAGE3(10,NPI) 

287 FORMAT (' ',F5.1,IX,10(2X,F5.1)) 

48 CONTINUE 

WRITE (IP, 297) CHAR(12) 

297 FORMAT (' ' ,A1) 

WRITE (IP, 307 ) HEAD,NP 
307 FORMAT (T6,'MAP OF DELTA VELOCITY AT LI 

* FOR ',A26,' TRAJECTORIES NODE OPTION ',11) 

WRITE (IP, 317 ) DATP,TIMP 

317 FORMAT (T27,'DATE',A10,' TIME',A10) 

WRITE (IP, 327 ) HEAD, FTIM, AINCEO * DPR, AINCM * DPR 




327 FORMAT (IX, A26,' FLIGHT TIME (HR) = ',F5.1,' INCL EARTH ', 

* F5.1,' INCL MOON = ',F5.1) 

WRITE (IP, 337) ALONO(1),ALONO(2),ALONO(3),ALONO(4),ALONO(5), 

*ALONO(6), ALONO(7), ALONO(8), ALONO(9), ALONO(10) 

337 FORMAT (/,' ALON> ', 10(2X, F5.1)) 

WRITE (IP,347) 

347 FORMAT (' ALAT') 

DO 58 NPI = 1 , 19 

WRITE (IP,357) ALATP(NPI),PAGE4(1,NPI),PAGE4(2,NPI), 

* PAGE4(3,NPI),PAGE4(4,NPI),PAGE4(5,NPI),PAGE4(6,NPI), 

* PAGE4(7,NPI),PAGE4(8,NPI),PAGE4(9,NPI),PAGE4(10,NPI) 

357 FORMAT(1X,F5.1,2X,10(1X,F6.0)) 

58 CONTINUE 

WRITE(IP,367) CHAR(12) 

367 FORMAT (' ',A1) 

WRITE(IP,377) HEAD, NP 

377 FORMAT(T2,'DELTA VELOCITY AT EARTH FOR ',A26,' TRAJECTORY 

* NODE OPTION ',11) 

WRITE(IP,387) DATP, TIMP 

387 FORMAT(T27,'DATE ', A10,' TIME ', A10) 

WRITE(IP,397) HEAD, FTIM, AINCEO * DPR, AINCM * DPR 
397 FORMAT(' ',A26,' FLIGHT TIME (HR) = ',F6.1,' INCL EARTH ', 

*F6.1,' INCL MOON = ',F6.1) 

WRITE(IP,407) ALONO(1), ALONO(2), ALONO(3), ALONO(4), 

*ALONO(5), ALONO(6), ALONO(7), ALONO(8), ALONO(9), ALONO(10) 

407 FORMAT(/,' ALON> ', 10(2X, F5.1)) 

WRITE(IP,417) 

417 FORMAT(' ALAT') 

DO 68 NPI = 1, 19 

WRITE(IP,427) ALATP(NPI), PAGE5(1,NPI), PAGE5(2,NPI), 

* PAGE5(3,NPI), PAGE5(4,NPI), PAGE5(5,NPI), PAGE5(6,NPI), 

* PAGE5(7,NPI), PAGE5(8,NPI), PAGE5(9,NPI), PAGE5(10,NPI) 

427 FORMAT (IX,F5.1,2X,10(IX,F6.0)) 

68 CONTINUE 

WRITE(IP,437) CHAR(12) 

437 FORMAT (' ' ,A1) 

WRITE(IP,447) HEAD, NP 

447 FORMAT(T2,'MAP OF AZMM FOR ',A26,'TRAJECTORY 

* NODE OPTION ',11) 

WRITE(IP,457) DATP, TIMP 

457 FORMAT(T27,'DATE',A10,' TIME',A10) 

WRITE(IP,467) HEAD, FTIM, AINCEO * DPR, AINCM * DPR 
467 FORMAT(IX,A26,' FLIGHT TIME (HR) = ',F5.1,' INCL EARTH ', 

*F5.1,' INCL MOON = ',F5.1) 

WRITE(IP,477) ALONO(1), ALONO(2), ALONO(3), ALONO(4), ALONO(5), 

* ALONO(6), ALONO(7), ALONO(8), ALONO(9), ALONO (10) 

477 FORMAT(/,' ALON> ', 10(2X, F5.1)) 

WRITE(IP,487) 

487 FORMAT (' ALAT') 




DO 78 NPI = 1, 19 

WRITE(IP,497) ALATP(NPI), PAGE6(1,NPI), PAGE6(2,NPI), 

* PAGE 6(3,NPI), PAGE 6(4,NPI), PAGE 6(5,NPI), PAGE6(6,NPI), 

* PAGE 6(7,NPI), PAGE 6(8,NPI), PAGE6(9,NPI), PAGE 6(10,NPI) 


497 

FORMAT(IX,F5. 

1,2X,10(IX,F6.0)) 

78 

CONTINUE 
GOTO 3000 
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CONTINUE 

COSALAT = 

QCOS 

(ALAT ) 


SINALAT = 

QSIN 

(ALAT ) 


COSALON = 

QCOS 

(ALON ) 


SINALON = 

QSIN 

(ALON ) 


COSPHI = 

COSALAT * COSALON 


COSAINC = 

QCOS 

(AINCM / DPR ) 


SINAINC = 

QSIN 

(AINCM / DPR ) 


RMR = COSPHI +QSQRT (COSPHI ** 2. - (1. - CMUE / CMUM )) 

RRM = RREM * 0.15 

C CALCULATION OF XYZ COORDINATES RELATIVE TO EARTH AT SPHERE 

C OF INFLUENCE 

C SPHERE OF INFLUENCE HAS BEEN REDEFINED EQUAL TO THE LIBRATION 

C POINT #1 RADIUS 

XXM = RRM * COSALAT * COSALON 

XX = -XXM + RREM 

YY = -RRM * COSALAT * SINALON 

ZZ = RRM * SINALAT 

YYM = -YY 

C WRITE (IS,507) XX,YY,ZZ 

C 507 FORMAT (' XYZ POSITION AT SPHERE ',Fll.0,IX,Fll.0,IX,Fll.0) 
RRE = QSQRT(XX**2.+YY**2.+ZZ**2. ) 

COSANGA = XX / RRE 

SINANGA = QSQRT (YY ** 2. + ZZ ** 2. ) / RRE 
C ANGA = QATAN (SINANGA / COSANGA ) 

ALONX - QATAN (YY / XX ) 

ALATX = QATAN (ZZ / QSQRT (XX ** 2. + YY **2. ) ) 

ANODEE = ALONX - QASIN (QTAN (-ALATX ) / QTAN (AINCE ) ) 

C END SPHERE OF INFLUENCE CALC 

30 CONTINUE 

ANODEM = ALON-AVM 

IF(ANODEM .LT.0.) ANODEM = ANODEM + 2.*PI 

CALL VELTRANS(VELM,GAMAM,AZMM,ALAT,ALON,VXM,VYM,VZM,VMAGM,DPR) 

VXM = VXM+XDLO 

VYM = VYM+YDLO 

QQEMIN = 2.1*RPE/(RRE+RPE) 

VELEMIN = QSQRT(QQEMIN*CMUE/RRE) 

IF(VELE .LT. VELEMIN) VELE = VELEMIN 
TIEMS = TIEM 

C TIMEES CHANGED TO TIEMS, TIMEE CHANGED TO TIEM 

35 CONTINUE 

CALL GAMACALC(RPE,VELE,RRE,CMUE,COSGAME,VPE,VCIRE, 


*TIME1, TRAJE, DPR, ALAT, ALON, FTNM) 

VELE2 = VELE + 10. 

CALL GAMACALC(RPE,VELE2,RRE,CMUE,COSGAME,VPE,VCIRE, 

*TIME2, TRAJE, DPR, ALAT, ALON, FTNM) 

DELVELE = 10./(TIME2-TIME1)*(FTIM-TIME1-TIMM) 

IF(DELVELE .LT. 500.) THEN 
VELE =VELE + DELVELE 
ELSE 

VELE = VELE + DELVELE/QABS(DELVELE)*500. 

END IF 

CALL GAMACALC(RPE,VELE,RRE,CMUE,COSGAME,VPE,VCIRE,TIEM, 

*TRAJE,DPR,ALAT,ALON,FTNM) 

IF(TIEM .EQ. 0.) GOTO 60 
TIMT = TIEM+TIMM 

IF(ABS(FTIM-TIMT) .GT. 1.) GOTO 35 

GAMAE = QFLOAT(MOOD)*QATAN(QSQRT(1.-COSGAME**2.)/COSGAME) 

ALONE = 180. / DPR + ALONX 

ALATE = QATAN(ZZ/QSQRT(XX**2.+YY**2.)) 

50 CONTINUE 

AINCE = AINCEO 

C 'CALCULATION OF EARTH ORBIT AZIMUTH AT SPHERE OF INFLUENCE 

IF(AINCE .NE. 0.0) THEN 

PHIE = PI-QASIN(ZZ/(RRE*QSIN(AINCE))) 

AZME = QATAN2(1./QTAN(AINCE),QCOS(PHIE)) 

ELSE 

AZME - PI/2 
END IF 
C 

CALL VELTRANS(VELE,GAMAE,AZME,ALATE,ALONE,VXE,VYE,VZE,VELE, 
* DPR) 

DELVEL = QSQRT ( (VXE-VXM) **2 . + (VYE-VYM) **2 . + (VZE-VZM) **2 . ) 
DVCIRE = VPE - VCIRE 

DVE = QSQRT (VCIRE ** 2. + VPE**2.-2. * VCIRE*VPE* 

*COS (AINCEO - AINCE ) ) 

DVTSAV = DVTOTAL 

DVTOTAL = DELVEL + DELVLP1 + DVCIRE 
C CALCULATE POSITION AND VELOCITY OF MOON 

C ITTERATE FOR FLIGHT TIM TO SPHERE OF INFLUENCE 

TIM = (TIMJ-2451545.)/36525.+TIEM/876600. 

CALL POSVELMO(TIM,RREMER,XDLO,YDLO) 

RREM = RREMER*20295741. 

IF(ABS(TIEM-TIEMS) .GT. .1) GOTO 25 

IF(DVTOTAL.LT. DVMIN)THEN 
DVMIN = DVTOTAL 
ALONMIN = ALON 
ALATMIN = ALAT 
VELMMIN=VELM 
VXMP = VXM 
VYMP = VYM 




VZMP = VZM 
VMAGMP = VMAGM 
GAMAMP = GAMAM 
AZMMP = AZMM 
AINCMP = AINCM 
ANODEMP = ANODEM 
TIMMP = TIMM 
DELVLP1P = DELVLP1 
DELVELP = DELVEL 
VXEP = VXE 
VYEP = VYE 
VZEP = VZE 
VELEP = VELE 
GAMAEP = GAMAE 
AZMEP = AZME 
AINCEP = AINCE 
ANODEEP = ANODEE 
TIEMP = TIEM 
DVCIREP = DVCIRE 
DVTOTALP = DVTOTAL 
END IF 

IF(DELVEL .LT. DVSIM)THEN 
DVSIM = DELVEL 
ALONSIM = ALON 
ALATSIM = ALAT 
VELMSIM = VELM 
OMEGE = PHIE - THETAE 
END IF 

60 CONTINUE 

WRITE (IS,517) ALON*DPR, ALAT*DPR 
517 FORMAT (' ALON = ',F6.1,' ALAT = ',F6.1) 

IF (ICALL.EQ.l) GOTO 1000 
IF (ICALL.EQ.2) GOTO 2000 
3000 STOP 
END 

SUBROUTINE FLYBY (ALON, ALAT, RREM, RPM, GAMAM, VELM, VPM, AIM, 
*DELVLP1,AZMM,AVM,TIEMM,MOOD) 

IMPLICIT REAL * 16 (A-H,0-Z) 

IS = 5 

DPR = 57.29578 
PI = 3.141593 
CMUE = 1.407647E+16 
CMUM = 1.731400E+14 
FTNM = 6076.115 
REE = 20925741. 

REMO = 5.7039E+6 
OMEGM ■ .9582117947E-2 
FTM = .3048 



R1 = .15 * RREM 

PRINT *,' TIMEM INCL THETAM GAMAM VELM QQ EEM' 

CONTINUE 

LOCATION OF LP#1 IN LUNAR RELATIVE COORDINATES. LIBRATION 
POINT ROTATES WITH THE MOON AND THEREFORE IS EFFECTED BY 
FLYBY TIME 

XI = -R1 * QCOS(OMEGM*TIEMM) 

Y1 = -R1 * QSIN(OMEGM*TIEMM) * FLOAT(MOOD) 

Z1 = 0. 

R1 IS LP#1, R2 IS SPACECRAFT AT SPHERE OF INFLUENCE. SPHERE 
OF INFLUENCE IS DEFINED EQUAL TO RADIUS OF LP#1. THIS IS 
SLIGHTLY LARGER THAN THE ONE USED IN THE PLANE CHANGE PROGRAM. 
HOWEVER, THIS CHANGE SIMPLIFIES THE CALCULATION OF THE FLYBY 
TRAJECTORY 

X2 = R1 * QCOS(ALAT) * QCOS(ALON) 

Y2 = R1 * QCOS(ALAT) * QSIN(ALON) 

Z2 = R1 * QSIN(ALAT) 

CALCULATION OF TRUE ANOMALY (THETA) 

DOT = (X1*X2 + Y1*Y2 + Z1*Z2) 

COS2THETA = DOT/R1/R1 

THETA = PI - QACOS(COS2THETA)/2. 

COSTHETAM = QCOS(THETA) 

CALCULATION OF ECCENTRICITY OF FLYBY ORBIT FROM RP, Rl, AND 
THETA 

EEM = (RPM/R1 - 1.)/(COSTHETAM - RPM/R1) 

COSGAMA = (1.+EEM*COSTHETAM)/QSQRT(1.+2.*EEM*COSTHETAM+EEM*EEM) 
CALCULATE FLIGHT PATH ANGLE. GAMA IS POSITIVE AT LI AND NEGA¬ 
TIVE AT SOI FOR OUTBOUND AND REVERSED FOR INBOUND TRAJECTORIES 

GAMAM = QACOS(COSGAMA) * FLOAT(MOOD) 

PPM = RPM * (1. + EEM) 

VPM = QSQRT(CMUM * (2./RPM + (EEM*EEM - 1.)/PPM)) 

VELM = QSQRT(CMUM * (2./Rl + (EEM*EEM - 1.)/PPM)) 

Rl CROSS R2 TO CALCULATE INCLINATION FROM Z3 
X3 = Y1*Z2 - Z1*Y2 
Y3 = Z1*X2 - X1*Z2 
Z3 = X1*Y2 - Y1*X2 

AMAG = QSQRT(X3*X3 + Y3*Y3 + Z3*Z3) 

AIM = QACOS(Z3/AMAG) 

QQM = Rl * VELM*VELM/CMUM 

COSAEX = (EEM + COSTHETAM)/(1. + EEM * COSTHETAM) 

AAX = Rl/(2. - QQM) 

TIEMMS = TIEMM 

IF (QQM .GT. 2.) GOTO 101 

CALCULATE FLIGHT TIME FOR ELLIPTIC ORBITS 
IF (ERRRP .GT. 0.) THEN 
TIMX = 0. 

GOTO 199 
END IF 




AEX = QACOS(COSAEX) 

SINAEX = QSQRT(1. - COSAEX*COSAEX) 

TIMX = AAX*QSQRT(AAX/CMUM)*(AEX-EEM*SINAEX)/3600. 

GOTO 199 
101 CONTINUE 

C CALCULATE FLIGHT TIME FOR HYPERBOLIC ORBITS 

COSHF = COSAEX 

SINHF = QSQRT(COSHF*COSHF - 1.) 

FFX = QLOG(COSHF + QSQRT (COSHF*COSHF - 1.)) 

TIMX = -AAX*QSQRT(-AAX/CMUM)*(EEM*SINHF - FFX)/3600. 

199 CONTINUE 

TIEMM = TIMX * 2. 

C CALCULATION OF AZIMUTH AND ANGLE AV AT R2 

ANODE = OMEGM * TIEMM * FLOAT(MOOD) 

IF (ALAT*FLOAT(MOOD) .GE. 0.) ANODE = ANODE + PI 
AVM = ALON - ANODE 

COSPHIM = QCOS(ALAT) * (QCOS(ANODE) * QCOS(ALON) + 

* QSIN(ANODE) * QSIN(ALON)) 

APHIM = QACOS(COSPHIM) 

IF (QABS(ALAT) .GT. 89.9/DPR) THEN 

AZMM = (PI + PI * ALAT/QABS (ALAT)) /2 . 

GOTO 200 
END IF 

SINAZM = QCOS(AIM)/QCOS(ALAT) 

COSAZM = QSIN(AIM) * QCOS(APHIM)/QCOS(ALAT) 

AZMM = QATAN2(SINAZM, COSAZM) 

200 CONTINUE 

C CALCULATE AZIMUTH FOR LI 

ALATL1 = 0. 

ALONL1 = QATAN2(Yl, XI) 

COSPHIL1 = QCOS(ALATLl) * (QCOS (ANODE) * QCOS(ALONLl) + 

* QSIN(ANODE) * QSIN(ALONL1)) 

APHIL1 = QACOS(COSPHIL1) 

SINAZM = QCOS(AIM)/QCOS(ALATLl) 

COSAZM = QSIN(AIM) * QCOS(APHIL1)/QCOS(ALATLl) 

AZMLP1 = QATAN2(SINAZM, COSAZM) 

C PRINT DATA AND ITERATE FOR FLIGHT TIME 

IF (QABS(TIEMM - TIEMMS) .GT. 0.1) GOTO 10 
C WRITE (IS,17) TIEMM, AIM*DPR, THETA*DPR, GAMAM*DPR, 

C *VELM, QQM, EEM 

C 17 FORMAT (IX,F5.1,3(2X,F5.1),2X,F10.2,2X,F10.4,2X,F10.6) 
IF (QABS(TIEMM - TIEMMS) .GT. 0.1) GOTO 10 
C LOCATION OF LB#1 IN LUNAR RELATIVE COORDINATES 

ALATLP1 = 0. 

ALONLP1 = 180. 

CALL VELTRANS(VELM,GAMAM,AZMLP1,ALATLP1,ALONLP1,VXLP1, 

* VYLP1,VZLP1,VMAGLP1) 

VYCOR = -OMEGM*Rl/3600. 

DELVLP1 = QSQRT(VXLP1**2.+(VYLP1 - VYCOR)**2.+VZLP1**2.) 




C GAMA IS NEGATIVE AT SOI, POSITIVE AT LI FOR OUTBOUND TRA- 

C JECTORIES, REVERSED FOR INBOUND FLIGHT. NEED TO CHANGE 

C SIGN FOR SOI CALCULATIONS 

GAMAM = -GAMAM 
RETURN 
END 

SUBROUTINE XYZPOS(RRX,ALAT,ALON,XX,YX,ZX) 

IMPLICIT REAL * 16 (A-Z) 

XX = -RRX*COS(ALAT)*COS(ALON)+RREM 
YX = -RRX*COS(ALAT)*SIN(ALON) 

ZX = RRX*SIN(ALAT) 

RETURN 

END 


SUBROUTINE POSVELMO(TIM,RRM,XDLO,YDLO) 
IMPLICIT REAL * 16 (A-H,0-Z) 

DELT = 0.5/36525./24./3600. 

T1 = TIM-DELT 

CALL MOON (Tl, RAM, DECM,RM1) 

T2 = TIM+DELT 

CALL MOON(T2,RAM,DECM,RM2) 

XDLO = (RM2-RM1)*20925741. 

RRM = (RM2+RM1)/2. 

RRM = RRM 

RRMB = RRM -7.412789E-01 
YDLO = 200570.2/RRMB 
RETURN 
END 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 


SUBROUTINE MOON (T, RAM, DECM, RM) 

FINDS LOCATION OF MOON IN EQUATORIAL COORDS. AT ANY TIM 
REF: '87 ASTRONOMICAL ALMANAC 

T IS JULIAN CENTURIES SINCE YEAR 2000 

LAM IS MOON'S ECLIPTIC LONGITUDE 

BETA IS MOON'S ECLIPTIC LATITUDE 

PIE IS HORIZONTAL PARALLAX 

RM IS DIST. TO MOON IN EARTH RADII 

RAM IS RT. ASCENSION OF MOON 

DECM IS MOON'S DECLINATION 

SD IS SEMIDIAMETER OF MOON'S ORBIT 


IMPLICIT REAL * 16 (A-Z) 
C PRINT *, ' MOON' 




P = 3.1415926535 
C = P / 180 

LAM = C*218.32+C*481267.883*T+C* 6.29 * QSIN(C * 134.9 + C 
*477198.85 * T) - C * 1.27 * QSIN(C * 259.2 - C * 413335.38 
*T) + c * .66 * QSIN(C * 235.7 + c * 890534.23*T) 

LAM = LAM + c * .21 * QSIN(c * 269.9 + c * 954397.7*T) - c 
*.19 * QSIN(c * 357.5 + c * 35999.05 * T) - c * .11 * 
*SIN(c * 186.6 + c * 966404.05*T) 

beta = c*5.13*QSIN(c*93.3 + c * 483202.03 * T) + c * 

* .28 * QSIN(C * 228.2 + C * 960400.87 * T) -c*.28*QSIN 

* (c*318.3+c* 60003.18*T)-c*.17*QSIN(c*217.6-c*407332.2 * T) 

pie = c*.9508+c*.0518*COS(c*134.9 + c * 477198.85 * T) + 
*c*.0095*QCOS(c*259.2-c*413335.38*T)+c* .0078*COS(c * 23 
*5.7+c*890534.23*T)+c*.0028*QCOS(c*269.9+c*954397.7* T) 

SD = .2725 * pie 
RM = 1. / QSIN(pie) 

1 = QCOS(beta) * QCOS(LAM) 

M = .9175 * QCOS(beta) * QSIN(LAM) - .3978 * QSIN(beta) 

n = .3978 * QCOS(beta) * QSIN(LAM) + .9175 * QSIN(beta) 

RAM = QATAN2(M, 1) 

DECM = QASIN(n) 

RETURN 

END 

SUBROUTINE GAMACALC (RPX, W, RRX, CMUX, COSGAMX, VPX, VCIRX, 
*TIMX, TRAJ, DPR, ALAT, ALON, FTNM) 

IMPLICIT REAL * 16 (A-H, O-Z) 

CHARACTER*10 TRAJ 

IHYPER = 1 

TRAJ = 'ELIPT ' 

QQX = RRX*W**2/CMUX 

IF(QQX-2 .LT. 1.0E-06) QQX = QQX - 1.0E-06 
IF (QQX .GT. 2.) THEN 
IHYPER = -1 
TRAJ = 'HYPER ' 

PRINT *,' ***** TRAJECTORY IS HYPERBOLIC ' 

END IF 

AAX = RRX/(2.-QQX) 

IF (AAX. GT . 1.0E12 .OR. AAX. LT. -1.0E12) AAX = -1.0E12 

EEX = 1.-RPX/AAX 

PPX = AAX*(1. -EEX ** 2.) 

COSGAMX = QSQRT(RPX/RRX*(1+EEX)/QQX) 

IF (COSGAMX .GT. 1.) THEN 



TIMX = 0. 

RETURN 
END IF 

GAMAX = QACOS(COSGAMX) 

VPX = QSQRT(CMUX*(1+EEX)/RPX) 

VCIRX = QSQRT(CMUX/RPX) 

COSTHETAX = (PPX/RRX-1)/EEX 

THETAX = QACOS(COSTHETAX) 

COSAEX = (EEX+COSTHETAX)/(l+EEX*COSTHETAX) 

ERRRP = RRX-AAX*(1+EEX) 

IF(ERRRP .GT. 0. .AND. AAX .GT. 0.) THEN 

WRITE (IS,537) ERRRP/6076.1155, ALAT*DPR, ALON*DPR 
537 FORMAT (' RADIUS > APOGEE BY NMI ',F7.5, F8.0, F7.5) 
WRITE (IS,547) QQX, AAX/FNTM, EEX, GAMAX*DPR, VPX 
547 FORMAT (' QQX = ',F7.5,' AAX = ',F8.0,' EEX = ',F7.5, 

* ' GAMAX = ', F5 .If' VPX = ',F7.1) 

END IF 

IF(QQX .GT. 2.) GOTO 101 
C CALC FLIGHT TIM FOR ELIPTICAL ORBITS 

IF (ERRRP .GT. 0.) THEN 
TIMX = 0. 

GOTO 199 
END IF 

AEX = QACOS(COSAEX) 

SINAEX = QSQRT(1-COSAEX**2.) 

TIMX = QSQRT(AAX ** 3. / CMUX)*(AEX-EEX*SINAEX)/3600. 
GOTO 199 
101 CONTINUE 

C CALC FLIGHT TIM FOR HYPERBOLIC ORBITS 

COSHF = COSAEX 
SINHF = QSQRT(COSHF**2.-1.) 

FFX = QLOG(COSHF+QSQRT(COSHF**2.-1.)) 

TIMX = QSQRT(-1. *AAX*AAX*AAX/CMUX)*(EEX*SINHF-FFX)/3600. 
199 CONTINUE 
RETURN 
END 


SUBROUTINE VELTRANS(VEL,GAMA,AZM,ALAT,ALON,VXX,VYX,VZX,VMAG 
k DPR) 

IMPLICIT REAL * 16 (A-Z) 

RRD = VEL * QSIN (GAMA ) 

RPHID = VEL * QCOS (GAMA ) 

VLON = RPHID * QSIN (AZM ) 

VLAT = RPHID * QCOS (AZM ) 

VXR = -RRD * QCOS (ALAT ) * QCOS (ALON ) 

VYR = -RRD * QCOS (ALAT ) * QSIN (ALON ) 

VZR = RRD * QSIN (ALAT ) 




VXLA = VLAT * QSIN (ALAT ) * QCOS (ALON ) 

VYLA = VLAT * QSIN (ALAT ) * QSIN (ALON ) 

VZLA = VLAT * QCOS (ALAT ) 

VXLO = VLON * QSIN (ALON ) 

VYLO = -VLON * QCOS (ALON ) 

VZLO =0.0 

VXX = VXR + VXLA + VXLO 
VYX = VYR + VYLA + VYLO 
VZX = VZR + VZLA + VZLO 

VMAG = QSQRT (VXX **2. + VYX ** 2.+ VZX ** 2. 

RETURN 

END 


Appendix C. Program Variables 


INPUT 

VARIABLE DESCRIPTION _ 

AINCEO Earth orbit inclination (degrees). This is the angle between the plane of the low 
Earth orbit and the plane of the Moon’s orbit about the Earth. 

ALATI Initial Sphere of Influence longitude for map (degrees) 

ALONI Initial Sphere of Influence latitude for map (degrees) 

DELLAT Incremental latitude for map (degrees) 

DELLON Incremental longitude for map (degrees) 

FITM Flight time for trajectory (hours) 

HPE Holding perigee altitude of Earth orbit (nautical miles) 

HPM Perigee altitude of Lunar "flyby" trajectory (nautical miles) 

MD Leg of trip for which to perform calculations (OUTBOUND or RETURN) 

TIMJ Earth departure Julian date (where January 1, 2000 is day 2,451,545. Refer to 

Section C of "The Astronomical Almanac of the Year 1988"). 

CONSTANT VALUE DESCRIPTION _ 


c 

n/180 

Degrees per radian (deg./rad.) 

CMUE 

1.407647E+16 

Gravitational parameter of the Earth (ftVsec 2 ) 

CMUM 

1.731432E+14 

Gravitational parameter of the Moon (ft 3 /sec 2 ) 

DPR 

57.29578 

Degrees per radian (deg./rad.) 

FTM 

0.3048 

Meters per foot (m/ft) 

FTNM 

6,076.115 

Feet per nautical mile (ft/nmi) 

P 

3.1415927 

n (dimensionless) 

PI 

3.141593 

11 (dimensionless) 

REE 

20,925,741 

Radius of the Earth (ft) 

REMO 

5,703,900 

Radius of the Moon (ft) 

RREM 

1,261,152,353 

Earth-Moon distance of centers (ft) 



AAX 


Semi-major axis of one of the transfer orbits 
AEX Eccentric anomaly of one of the transfer orbits 

AIM Inclination of "flyby" trajectory 

AINC LI orbit inclination in radians (always 0) 

AINCE Earth orbit inclination in radians 

AINCEP Earth orbit inclination in radians 

AINCM "Flyby" orbit inclination (degrees) 

AINCMP "Flyby" orbit inclination (degrees) for the SOI point associated with minimum 
AV. 

ALAT Initial latitude for map (radians) 

ALATE Latitude of SOI point, measured from the Earth 

ALATL1 Latitude of LI, measured from the Moon (0°) (Moon-fixed) 

ALATLP1 Latitude of LI, measured from the Moon (0°) 

ALATMIN Latitude of SOI point associated with the minimum total AV 
ALATPO ALAT in degrees 

ALATSIM Latitude of SOI point associated with the minimum SOI AV 
ALATX Latitude of the SOI point, measured from the Earth 
ALON Initial longitude for map (radians) 

ALONE Longitude of SOI point, measured from zero longitude at Earth (-X direction) 

ALONL1 Longitude of LI, measured from the Moon (Moon-fixed) 

ALONLP1 Longitude of L1, measured from the Moon (180°) 

ALONMIN Longitude of SOI point associated with the minimum total AV 
ALONO() Variable that contains the print matrix column headings (longitudes) 

ALONSIM Longitude of SOI point associated with the minimum SOI AV 


ALONX 


AMAG 

ANGA 

ANODE 

ANODEE 

ANODEEP 

ANODEM 

ANODEMP 

APHIL1 

APHIM 

AVM 

AZM 

AZME 

AZMEP 

AZMLP1 

AZMM 

AZMMP 

BETA 

CMUX 

COSAEX 

COSAINC 

COSALAT 


Longitude of the SOI point, measured from the Earth-Moon line at the Earth 

Magnitude of position vector cross-product (R2 X Rl) used to calculate "flyby" 
inclination 

The angle between the Earth-to-Moon line and the Earth-to-SOI-point line 
The longitude of the ascending or decending node of "flyby" trajectory 
Longitude of the Earth ascending or descending node 

Longitude of the Earth ascending or descending node associated with minimum 
AV 

Longitude of the Lunar ascending or descending node 

Longitude of the Lunar ascending or descending node associated with the 
minimum AV 

Angle between line of nodes and LI 

Angle between line of nodes and Moon-to-SOI line 

Angle between the Lunar node and the projection of the SOI point onto the Earth- 
Moon plane (D) 

Azimuth of the SOI point 

Azimuth (from the Earth) of the SOI point for Earth-SOI trajectory 

Azimuth (from the Earth) of the SOI point for Earth-SOI trajectory associated 
with minimum AV 

Azimuth of "flyby" orbit at LI 

Azimuth (from the Moon) of the SOI point for "flyby" trajectory 

Azimuth (from the Moon) of the SOI point for "flyby" trajectory with minimum 
AV 

Moon’s ecliptic latitude 

Earth or Moon gravitational parameter 

Cosine of the eccentric anomaly of one of the transfer orbits 

Cosine of the LI orbit inclination (always 1) 

Cosine of the SOI point latitude 




COSALON Cosine of the SOI point longitude 

COSANG A Cosine of the angle between the Earth-to-Moon line and the Earth-to-SOI-point 
line 

COSAZM Cosine of the azimuth angle at the SOI point or LI for the "flyby" trajectory 

COSGAMA Cosine of flight path angle at SOI for "flyby" trajectory 

COSGAME Cosine of the flight path angle at SOI of the Earth-to-SOI trajectory 

COSGAMX Cosine of the flight path angle at SOI of one of the transfer orbits 

COSHF Hyperbolic cosine of the eccentric anomaly of one of the transfer orbits 

COSPHI Cosine of the angle between the Moon-to-Earth line and the Moon-to-SOI-point 
line 

COSPHIL1 Cosine of angle between line of nodes and LI 

COSPHIM Cosine of the angle between the line of nodes and the Moon-to-SOI line 

COSTHETAM Cosine of the true anomaly of the "flyby" orbit 
COSTHETAXCosine of the true anomaly of one of the transfer orbits at SOI 
COS2THETA Cosine of two times the true amomaly of the "flyby" 

DATP Today’s date 

DECM Declination of the Moon 

DELT A fraction of TIM that represents a half-second 

DELV(,) Total AV for matrix of SOI longitudes and latitudes 

DELVEL AV at the SOI point to patch the "flyby" trajectory with the Earth-to-SOI 
trajectory 

DELVELP AV at the SOI point to patch the "flyby" trajectory with the Earth-to-SOI 
trajectory associated with minimum AV 

DELVELE Estimated AV required for the Earth-to-SOI trajectory in order to satisfy the time- 
of-flight requirement 

DELVLP1 AV for circularization at LI 

DELVLP1P AV for circularization at LI associated with minimum AV 



DOT 

DVCIRE 

DVCIREP 

DYE 

DVMIN 

DVSIM 

DVTOTAL 

DVTOTALP 

DVTSAV 

EEM 

EEX 

ERRRP 

FFX 

GAMA 

GAMAE 

GAMAEP 

GAMAM 

GAMAMP 

GAMAX 

HEAD 

ICALL 

II 

IHYPER 

IP 


Dot product of position vectors for LI and SOI point (R1 R2) 

AV between Earth circular orbit and Earth-to-SOI trajectory 

AV between Earth circular orbit and Earth-to-SOI trajectory associated with 
minimum AV 

Unused plane change variable 
Hold variable for minimum total AV 
Hold variable for minimum SOI AV 
Total AV for flight 

Total AV for flight associated with minimum AV 

Temporary storage for total AV 

Eccentricity of the "flyby" orbit 

Eccentricity of one of the transfer orbits 

Difference between orbital range at SOI and apogee range 

Hyperbolic eccentric anomaly 

Flight path at SOI of one of the transfer orbits 

Flight path at SOI of Earth-to SOI trajectory 

Flight path at SOI of Earth-to SOI trajectory associated with minimum AV 

Flight path at SOI of "flyby" trajectory 

Flight path at SOI of "flyby" trajectory with minimum AV 

Flight path angle at SOI of one of the transfer orbits 

Heading that indicates leg of journey (e.g., "Earth to LI Flyby" for outbound) 

Flag to track occurance of call to in-program subroutine. 

Counter (1-19), representing increments in latitude for the output matrices 
Indicator that describes whether an orbit is hyperbolic 
Output number indicating output to file 




IPRINT 


Flag for print block that contains matrix column headers 
IS Output number indicating output to screen 

L A geocentric direction cosine 

LAM Moon’s ecliptic longitude 

M A geocentric direction cosine 

MOOD Leg of journey (+1: outbound) 

N A geocentric direction cosine 

NN Counter (1-10) representing increments in longitude for the output matrices 

NPI Counter for the 19 matrix rows during report printing 

OMEGE Argument of perigee for Earth-SOI trajectory 
OMEGM Angular velocity of Moon 

PAGE1() Total AV for a given SOI longitude and latitude. Cell values for Report #1 matrix 

PAGE2(,) SOI AV for a given SOI longitude and latitude. Cell values for Report #2 matrix 

PAGE3(,) Inclination of "flyby" trajectory for a given SOI longitude and latitude. Cell 

values for Report #3 matrix 

PHIE Used to determine the azimuth of the SOI point 

PIE Horizontal parallax 

PPM Semi-latus rectum of the "flyby" orbit 

PPX Semi-latus rectum of one of the transfer orbits 

QQEMIN Vis-viva parameter for the Earth-to-SOI-point trajectory 

QQM Vis-viva parameter for the "flyby" orbit 

QQX Vis-viva parameter for an orbit 

R1 Distance of LI and SOI surface from center of Moon (15% of Earth-Moon 

distance) 

RAM Right ascension of the Moon 

RM Distance from the Earth to the Moon in Earth radii 



RMR 

RM1 

RM2 

RPE 

RPHID 

RPM 

RPX 

RRD 

RRE 

RREMER 

RRM 

RRMB 

RRX 

SD 

SINAEX 

SINAINC 

SINALAT 

SINALON 

SINANGA 

SINAZM 

SINHF 

T 

THETA 

THETAX 

TIEM 


Ratio of Earth-Moon distance to Moon-SOI distance 
Distance from the Earth to the Moon in Earth radii 
Distance from the Earth to the Moon in Earth radii 
Distance from Earth’s center to Earth perigee orbit 
Orbital path component of the velocity vector 
Distance from Moon’s center to Lunar "flyby" perigee 
Distance from Earth center to perigee 
Radial component of the velocity vector 
Distance from the Earth to the SOI point 
Moon’s distance from the Earth (Earth radii) 

Distance from the Moon to the SOI point, when the Moon is at its mean distance 
from the Earth 

Distance of the Moon from the Earth-Moon baricenter 

Distance from Earth to SOI 

Semi-diameter of the Moon’s orbit 

Sine of the eccentric anomaly of one of the transfer orbits 

Sine of the Lunar orbit inclination 

Sine of the latitude of the SOI point 

Sine of the longitude of the SOI point 

Sine of the angle between the Earth-to-Moon line and the Earth-to-SOI point line 

Sine of the azimuth angle at the SOI point or LI for the "flyby" trajectory 

Hyperbolic sine of the eccentric anomaly of one of the transfer orbits 

Number of Julian centuries since the year 2000 AD 
True anomaly of "flyby" at the SOI 

True anomaly of one of the transfer orbits at SOI 

Earth-to-SOI time of flight (seconds) 




i i 11 i i p 


TIEMP 

TIEMS 

TIM 

TTME1 

TIME2 

TIEMM 

TIEMMP 

TIEMMS 


T2 

VCIRE 

VCIRX 

VEL 

VELE 

VELEP 

VELE2 

VELEMIN 

VELM 

VELMMIN 


Earth-to-SOI time of flight (seconds) associated with minimum AV 
Temporary storage for Earth-to-SOI time of flight 

Time of arrival at SOI point from Earth, in centuries since the year 2000 AD 
Earth-to-SOI time of flight (seconds) 

Earth-to-SOI time of flight (seconds) 

Time of flight from SOI to LI 

Time of flight from SOI to LI associated with minimum AV 

Temporary storage for SOI to LI time of flight 

Time Now 

Total time of flight 

Time of perigee passage for "flyby" 

Text that describes whether an orbit is hyperbolic or elliptical 

Text that describes whether the Earth-to-SOI trajectory is hyperbolic or elliptical 

Text that describes whether the "flyby" trajectory is hyperbolic or elliptical 

One-half second before TIM 

One-half second after TIM 

Velocity of Earth circular orbit 

Velocity of the Earth circular orbit 

Velocity at SOI of one of the transfer orbits 

Velocity at SOI of the Earth-to-SOI trajectory 

Velocity at SOI of the Earth-to-SOI trajectory associated with minimum AV 
Ten feet per second more than VELE 

Minimum velocity required such that apogee of the trajectory is just at SOI 
Velocity at SOI of the "flyby" trajectory 

Velocity at SOI of the "flyby" trajectory associated with the minimum total AV 




VELMOUT(, 

VELMSIM 

VLAT 

VLON 

VMAG 

VMAGLP1 

VMAGM 

VMAGMP 

VPE 

VPM 

VPX 

w 

VYMESTE 

VXE 

VXEP 

VXLA 

VXLO 

VXLP1 

VXM 

VXMP 

VXR 

VXX 

VYCOR 

VYE 


) Velocity at SOI of the "flyby" trajectory, for a given SOI longitude and latitude 
Velocity at SOI of the "flyby" trajectory associated with the minimum SOI AV 
Latitude component of the velocity vector 
Longitude component of the velocity vector 
Velocity vector magnitude 

Velocity vector magnitude at LI for "flyby" trajectory 

Velocity vector magnitude at SOI for "flyby" trajectory 

Velocity vector magnitude at SOI for "flyby" trajectory with minimum AV 

Perigee velocity of the Earth-to-SOI trajectory 

Perigee velocity of the "flyby" trajectory 

Perigee velocity for one of the transfer orbits 

Velocity at SOI of one of the transfer orbits 

Velocity at SOI point (apogee) of Earth-to-SOI trajectory 

X-coordinate of velocity vector at SOI for Earth-to-SOI trajectory 

X-coordinate of velocity vector at SOI for Earth-to-SOI trajectory associated with 
minimum AV 

X-component of the latitude component of the velocity vector 
X-component of the longitude component of the velocity vector 
X-component of the velocity vector at LI for "flyby" 

X-coordinate of velocity vector at SOI for "flyby" trajectory 

X-coordinate of velocity vector at SOI for "flyby" trajectory with minimum AV 

X-component of the radial component of the velocity vector 

Total X-component of the velocity vector 

Y-component of velocity for orbit at LI (no X- or Z- component) 

Y-coordinate of velocity vector at SOI for Earth-to-SOI trajectory 




VYEP 

VYLA 

VYLO 

VYLP1 

VYM 

VYMP 

VYR 

VYX 

VZE 

VZEP 

VZLA 

VZLO 

VZLP1 

VZM 

VZMP 

VZR 

VZX 

XI 

X2 

X3 

XDLO 

XX 

XXM 

Y1 


Y-coordinate of velocity vector at SOI for Earth-to-SOI trajectory associated with 
minimum AV 

Y-component of the latitude component of the velocity vector 

Y-component of the longitude component of the velocity vector 

Y-component of the velocity vector at LI for "flyby" 

Y-coordinate of velocity vector at SOI for "flyby" trajectory 

Y-coordinate of velocity vector at SOI for "flyby" trajectory with minimum AV 

Y-component of the radial component of the velocity vector 

Total Y-component of the velocity vector 

Z-coordinate of velocity vector at SOI for Earth-to-SOI trajectory 

Z-coordinate of velocity vector at SOI for Earth-to-SOI trajectory associated with 
minimum AV 

Z-component of the latitude component of the velocity vector 
Z-component of the longitude component of the velocity vector 
Z-component of the velocity vector at LI for "flyby" 

Z-coordinate of velocity vector at SOI for "flyby" trajectory 

Z-coordinate of velocity vector at SOI for "flyby" trajectory with minimum AV 

Z-component of the radial component of the velocity vector 

Total Z-component of the velocity vector 

X-component of the distance form the Moon to LI 

X-component of the distance from the Moon to the SOI penetration point 

X-component of cross-product (R1 X R2) 

X-coordinate of the velocity of the Moon 

Distance in the X-direction from the Earth to the SOI point’s X-coordinate 
Distance in the X-direction from the Moon to the SOI point’s X-coordinate 
Y-component of the distance form the Moon to LI 





Y2 Y-component of the distance from the Moon to the SOI penetration point 

Y3 Y-component of cross-product (R1 X R2) 

YDLO Y-coordinate of the velocity of the Moon 

YY Distance in the Y-direction from the Earth-Moon line to the SOI point’s 

coordinate 

YYM Distance in the Y-direction from the Moon to the SOI point’s Y-coordinate 

Z1 Z-component of the distance form the Moon to LI 

Z2 Z-component of the distance from the Moon to the SOI penetration point 

Z3 Z-component of cross-product (R1 X R2) 

ZZ Distance in the Z-direction from the Earth-Moon plane to the SOI point 




Appendix D. Detailed Program Description 

This section is a line-by-line description of the lines of code in the main program of LP1. Refer 
to Appendix B for the actual code listing, and figures K-l through K-3 for supporting illustra¬ 
tions. 

1. Declare the matrices DELV, VELMOUT, PAGE1, PAGE2, PAGE3, PAGE4, PAGE5, 
PAGE6, ALONO, and ALATP. 

2. Open the output file (LPl.OUT). 

3. Define the program constants. 

4. Read the program inputs. 

A. MD (leg of trip for which to perform calculations; outbound or return) 

B. HPE (perigee altitude of Earth circular orbit) 

C. HPM (perigee altitude of Lunar "flyby" trajectory) 

D. TIMJ (Earth departure Julian date. Default is 2451545, representing Jan. 1, 2000 
AD) 

E. ALONI, DELLON (Initial longitude and increment for output map) 

F. ALATI, DELLAT (Initial latitude and increment for output map) 

G. AINCEO (Angle between the plane of the low Earth orbit and the Earth-Moon 
plane) 

H. FTTM (flight time) 

5. Calculate the distance from the Earth’s center to Earth perigee orbit. 

RPE = HPE * FTNM + REE 

6. Calculate the distance from the Moon’s center to Lunar perigee orbit. 

RPM = HPM * FTNM + REMO 

7. Store the inclinations in temporary variables (A1NCE and AINC). 

8. Record the current date and time (DATP and TIMP). 

9. Print the headings for Report #1, "Velocity Map for Inbound/Outbound Trajectories". 

10. Echo back the input values onto Report #1. 

11. Initialize the velocity map matrix coordinates to 1,1 (II = rows, or latitude increments; 
NN = columns, or longitude increments). Set the flag IPRINT to zero, which has the 




effect of causing the longitude increments to be printed as a subheading. Initialize the 
minimum AV hold variable (DVMIN) to 99999. 

12. Calculate the orbital velocity of the Moon. 

YDLO = (iE/RREM 

13. Convert the initial map latitude and longitude from degrees to radians (ALAT and 
ALON). 

14. Initialize DVSIM at 99999. This variable will be used to hold the lowest AV for the 
trajectory, selected from all combinations of latitude and longitude. 

15. Initialize the current Total AV matrix cell (DELV) to 99999. 

16. Initialize with zeros the print line variables for the current cell on each report (PAGE1, 
PAGE2, PAGE3, PAGE4, PAGE5, PAGE6). 

17. Call the flyby subroutine to calculate the characteristics for a "flyby" trajectory between 
the current SOI point and LI. Pass input values for SOI latitude and longitude, Earth- 
Moon distance, and perigee altitude; and return output values for flight path angle and 
velocity (same for SOI point and LI), perigee velocity, time of flight, delta velocity 
needed to circularize at LI (outbound) or initiate the "flyby" (return), azimuth angle, 
inclination of the "flyby" trajectory, and the angle between the lunar nodes and the SOI 
point projection. 

18. Call the in-program subroutine given the current SOI conditions for the "flyby" trajectory 
(latitude, longitude, inclination, flight path, azimuth, and velocity) and iterate for a 
transfer orbit between LEO and the SOI point (with a velocity vector correction at the 
SOI) and determine the total AV and total flight time needed. 

19. Store the total AV, SOI AV, and "flyby" inclination in arrays for output. 

20. If all the matrix columns have been processed for this latitude, and if this has been the 
first row of the matrix, print the matrix column headings for Report #1. 

21. If all ten matrix columns have not been processed for this latitude, increment the column 
counter by one and increment the longitude by the amount specified by the user in the 
inputs. Return to step 15 to process the next longitude/latitude combination. 

22. If all ten matrix columns have been processed for this latitude, print the matrix row for 
this latitude for Report #1. 

23. If all 19 matrix rows have not been processed, increment the row counter by one, 
increment the latitude by the amount specified by the user in the inputs, reset the column 
counter to one, reset the longitude to the initial input longitude, and return to step 15 to 
process the next longitude/latitude combination. 

24. If all 19 matrix rows have been processed, continue with the following steps. 

25. Print the bottom section of Report #1, and the remaining five reports. 




Appendix E. In-Program Subroutine Description 


This section is a line-by-line description of the lines of code in the subroutine that is imbedded in 
program LP1, beginning at line #25. Refer to Appendix B for the actual code listing, and figures 
K-2 through K-4 for supporting illustrations. 

I. Establish the key distances and angles of the current SOI point from the Earth and Moon. 
A. Calculate the distance from the Moon to the SOI point. 

1. Define variables for the sine and cosine of the SOI latitude and longitude 
(COSALAT, COSALON, SINALAT, SINALON). 

2. Given a point on the SOI defined by the latitude and longitude, define the 
angle between the Moon-to-Earth line and the Moon-to-SOI line. 

cos( ) = cos(lat) * cos(lon) 

3. Define variables for the sine and cosine of the Lunar orbit inclination 
(COSAINC and SINAINC). 

4. Calculate the ratio of Earth-Moon distance to Moon-SOI distance 
(RMR = REM / RM). 

a. RRE 1 = RRM 2 + RREM 2 - 2 * RRM * RREM * cos( ) 

b. pE/RRE 2 = pM/RRM 2 .’. RRE 2 /RRM J = pE/pM 

(because p/R z = gravitational acceleration, and by definition, at the 
sphere of influence, acceleration towards the Earth equals 
acceleration towards the Moon). 

c. RRE7RRM 2 =pE/pM= 1+RREM 2 /RRM 2 -2*(RREM/RRM)*cos( ) 
= 1 + RMR 2 - 2*RMR*cos( ) 

d. 0 = -(pE/pM) + 1 + RMR 2 - 2*RMR*cos( ) 

= RMR 2 - 2*RMR*cos( ) +1 - (pE/pM) 

e. For the quadratic formula: 

a = 1, b = -2*cos( ), c = 1 - (pE/pM) 

RMR = (2*cos( )± 4cos 2 -4(1-pE/pM) )/2 
RMR = cos( ) + cos 2 -1 + pE/pM 

5. Determine the distance from the Moon to the SOI. 

RRM = RREM/RMR 




B. Determine an Earth-based rectangular coordinate system such that the X-direction 
is on the line from the Earth to the Moon; the Y-direction is in the direction of the 
Moon’s orbit; and the Z-direction is perpendicular to X and Y in right-hand 
coordinates. Project the SOI point onto the X-Y plane (Earth-Moon plane) in 
order to calculate the X- and Y-coordinates. 

1. Determine the distance from the Moon to the X-intercept. 

XXM = RRM * cos( ) = RRM * COSALAT * COSALON 

2. Determine the distance from the Earth to the X-intercept. 

XX = RREM - XXM 

3. Determine the Y-coordinate. Note that longitude increases in the negative 
Y direction. 

YY = -RRM * COSALAT * SINALON 

4. Determine the Z-coordinate. 

ZZ = RRM * SINALAT 

5. Determine the distance from the Moon to the Y-intercept. 

YYM = -YY 

C. Calculate the distance from the Earth to the SOI point. 

RRE = XX 2 + YY 2 + ZZ 2 

D. Identify the angle between the Earth-to-Moon line and the Earth-to-SOI line. 

1. cos(ANGA) = XX / RRE 

2. sin(ANGA) = YY 2 + ZZ 2 / RRE 

3. ANGA = tan'(sin(ANGA) / cos(ANGA)) 

E. Determine the Earth-based latitude and longitude of the SOI point. 

1. ALONX = tan'(YY / XX) 

2. ALATX = tan’ (ZZ/ XX 2 + YY 2 ) 

F. Determine the Earth ascending node for the orbit in which the trans-SOI bum will 
occur. 

1. Given the latitude of the bum at Earth (negative SOI latitude) and the 
inclination of the Earth orbit, spherical coordinate trigonometry states that 
the longitude from the node to the bum point is 
u = sin' 1 (tan(-ALATX) / tan(AINCE)) 


2. Longitude of the bum point (measured from the Earth-Moon line) is 
ALONX. 




3 . 


The node is calculated by subtracting the longitude in (a) above from the 
longitude in (b). 

ANODEE = ALONX - sin' 1 (tan(-ALATX) / tan(AINCE)) 

II. For the SOI-to-Ll trajectory, calculate Lunar node positions, and velocity vector at the 
current SOI point with respect to Earth. 

A. Determine the node. 

ANODEM = ALON - AVM 

B. Determine the components of the velocity vector at the SOI point. 

1. Call the subroutine VELTRANS to determine the X-, Y-, and Z- 
components of the velocity vector, from the viewpoint of the moon. 

2. Add the X- and Y-components of the Moon’s velocity to determine the 
SOI velocity from the viewpoint of the Earth. 

VXM = VXM + XDLO VYM = VYM + YDLO 

HI. For the Earth-to-SOI trajectory, calculate velocity at the SOI point, AV at Earth circular 
orbit, and time of flight. 

A. Modify the velocity at SOI, if necessary, to bring the trajectory up to at least the 
minimum velocity required to intercept the SOI. 

1. Calculate the vis-viva parameter QQEMIN for the Earth-to-SOI trajectory. 

QQEMIN = 2 - (RRE / a) where a = (RPE + RRE) / 2 
QQEMIN = 2 - ((2 * RRE) / (RPE + RRE)) 

QQEMIN = 2 * (1 - (RRE / (RPE + RRE))) 

QQEMIN = 2 * (RPE / (RPE + RRE)) 

2. Calculate the minimum velocity required, such that apogee is just at SOI. 

QQEMIN s (V 2 * RRE) / CMUE 

.-. VELEMIN = (QQEMIN * CMUE) / RRE 

3. Increase the velocity slightly so that the calculations converge. 

4. Compare the most recently calculated SOI velocity (VELE) with the 
minimum velocity. If it is too low, bring it up to the minimum. 

B. Temporarily store (until step IV-D) the calculated Earth-to-SOI time of flight. 

C. Further increase the velocity at SOI, if necessary, to meet the Earth-to-SOI time 
of flight requirement (but capping the AY at 500 ft/sec). Calculate the new SOI 
flight path angle, AV at Earth circular orbit, and time of flight. 

1. Call the subroutine GAMACALC to determine the time of flight from 
Earth perigee to SOI (TIMEE1), given the velocity VELE. 




2. Increase the velocity by 10 ft/sec, and call GAMACALC again to 
determine the new time of flight (TIMEE2). 

3. For the Earth-to-SOI leg, determine the ratio of time-of-flight shortfall to 
the time-of-flight increase that was due to the 10 ft/sec increase in 
velocity. 

Shortfall/Increase = (FTIME-TIMEE1 -TIMEM) / (TIMEE2 - 

T1MEE1) 

Apply this ratio to the increased velocity of 10 ft/sec to estimate the AV 
necessary to meet the required time of flight. 

DELVELE = 10 * (FTIME - TTMEE1 - TIMEM) / (TIMEE2 - 

TIMEEl) 

4. Cap this AV at 500 ft/sec, and add it to the original velocity. 

5. Call the subroutine GAMACALC to determine the new Earth-to-SOI time 
of flight for the revised velocity. 

D. Compare the total Eaith-to-Moon time of flight to the required time of flight. If 

the total is short, return to step M-C. 

E. Determine the components of the velocity vector at the SOI point. 

1. Calculate the flight path angle at SOI (GAMAE) for the Earth-to-SOI 
trajectory. The cosine of the flight path angle will be between zero and 
one, corresponding to angles between zero and 90. In reality, Moon-to- 
SOI flight path angles will range from zero to 90, but SOI-to-Moon angles 
will range from zero to -90. Gamma must be derived from cos(GAM) and 
multiplied by -MOOD (from the input: +1 for outbound and -1 for return) 
to get the correct flight path angle. The arctan function is used instead of 
arccos to evaluate gamma because of that function’s better debugging 
diagnostics capabilities. 

GAMAE = -MOOD * tan '( 1 - COSGAME 2 / COSGAME) 

2. Determine the longitude of the SOI from the Earth’s viewpoint. 

ALONE = ALONX +180 (since zero longitude is on the Earth- 

Moon line on the side of the Earth awav from the Moon). 

3. Determine the latitude of the SOI point. 

ALATE = tan' 1 (ZZ / XX 2 + YY 2 ) 

4. Calculate the azimuth of the SOI point 

AZME = tan' (l/tan(AINCE), cos(PHIE)) 

where PHIE = 11- sin '((ZZ/RRE)/sin(AINCE)) 

5. Call the subroutine VELTRANS to determine the X-, Y-, and Z-compo- 
nents of velocity at SOI for the Earth-to-SOI trajectory. 




IV. Combine the SOI-to-Moon calculations with the Earth-to-SOI calculations to identify 
total flight characteristics. 

A. Calculate total AV 

1. Calculate AV at the SOI point required to patch the Earth-to-SOI 
trajectory into the SOI-to-Ll trajectory. (This value is calculated by 
determining the difference between each of the components of the two 
trajectories). 

DELVEL = (VXE - VXM) 2 + (VYE - VYM) 2 + (VZE - VZM) 2 

2. Calculate AV at Earth perigee. 

DVCIRE = VPE - VCIRE 

3. Save the old total AV. 

DVTSAV = DVTOTAL 

4. Calculate the new total AV. 

DVTOTAL = DELVEL + DELVLP1 + DVCIRE 

B. Determine the distance from the Moon to the Earth for this trajectory just 
calculated. Use this value in subsequent iterations of the in-program subroutine, 
if necessary. 

1. Determine the time of arrival at SOI from Earth, in centuries since the 
year 2000 AD. 

TIM = Departure date + Time of flight - 2000 

where Departure date = TIMJ / 36525 days per century 
Time of flight = TIEM hrs / 876600 hrs per century 
The year 2000 = Day 2451545 / 36525 days per century. 

2. Call the subroutine POSVELMO to determine the Moon’s distance from 
the Earth and velocity at the time of arrival at SOI from the Earth. 

3. Convert Lunar distance from Earth radii to feet. 

RREM = RREMER * 20295741 ft/radii. 

C. Compare the value of Earth-to-SOI time of flight saved in step III-B above with 
the new value of time of flight calculated in step III-C above. If the difference 
has not yet converged (to less than 0.1 hour), reiterate this in-program subroutine 
beginning at step I-A. Otherwise, continue. 

D. If the newly calculated total AV is the smallest AV calculated so far, then store it 
in the variable DVMIN. Also store all the significant orbital characteristics of the 
trajectory for later use in the bottom section of Report #1. 

E. If the newly calculated SOI AV is the smallest calculated so far, then store it in 
the variable DVSIM. Also store its corresponding SOI longitude (ALONSIM), 
SOI latitude (ALATSIM), and SOI velocity for the SOI-to-Moon trajectory 
(VELMSIM). 




Appendix F. Subroutine FLYBY 


The subroutine FLYBY receives the parameters SOI longitude (ALON) and latitude (ALAT), 
Earth-Moon distance, "flyby" perigee radius (RPM), and flight mode (MOOD). FLYBY 
contains the same list of constants as in the main program, plus the radius of the SOI (Rl) and 
the mean angular speed of the Moon (OMEGM). It calculates the magnitude of the velocity at 
the SOI and at LI (VELM), the flight path angle at the SOI and at LI (GAMAM), the inclination 
of the "flyby" trajectory (AIM), the circularization AV needed at LI (DELVLP1), the azimuth 
angle at the SOI (AZMM), the angle between the Lunar node and the SOI projection (AVM), 
and the flight time between SOI and LI (TIEMM). 

As mentioned briefly in the introduction, the FLYBY routine takes advantage of one simplifi¬ 
cation. By enlarging the SOI radius to include LI from 11% of the Earth-Moon distance to 15%, 
the problem can be simplified considerably without significant error. Since the SOI penetration 
point is at the same distance from the Moon as LI, perigee passage will occur exactly half-way 
between the two points for any orbit containing them. Consequently, the true anomaly, flight 
path angle, and velocity will be the same at both the SOI and LI. 

1. Initialize constants. 

2. Determine location of LI in Lunar relative coordinates. Remember, LI rotates with the 
Moon, and therefore is affected by "flyby" time. X-axis is measured form Lunar center 
to Earth. Y-axis is measured in the opposite direction of the Moon’s velocity. 

XI = -Rl * cos(OMEGM * TIEMM) 

Y1 = -Rl * sin(OMEGM * TIEMM) * MOOD 
Z1=0 

3. Determine location of SOI point. 

X2 = Rl * cos(ALAT) * cos(ALON) 

Y2 = Rl * cos(ALAT) * sin(ALON) 

Z2 = Rl * sin(ALAT) 

4. Calculate true anomaly. 

A. Dot product of position vectors for SOI point and LI. 

DOT = XI * X2 + Y1 * Y2 + Z1 * Z2 

B. From geometry: 

cos(20) = DOT/R1/R1 

C. It follows: 

0 = arccos(cos(20))/2 

5. Calculate eccentricity of "flyby" orbit using general conic equation. 

EEM = (RPM/R1 - l)/(cos(0) - RPM/R1) 

6. Calculate flight path angle. 

cos(y) = (1 + EEM * cos(0))/ (1 + 2 * EEM * cos(0) + EEM 2 ) 
y= arcos(cos(Y)) * MOOD 




7. Calculate semi-latus rectum. 

PPM = RPM * (1 + EEM) 

8. Calculate perigee velocity. 

VPM = (m, * (2/RPM + (EEM 2 - 1)/PPM)) 

9. Calculate velocity at SOI and LI. 

VELM = (p m * (2/R1 + (EEM 2 - 1)/PPM)) 

10. Calculate inclination by crossing position components (R2 X R1). 

X3 = Y2 * Z1 - Z2 * Y1 
Y3 = Z2 * XI - XI * X2 
Z3 = X2 * Y1 - XI * Y2 
MAG = (X3 * X3 + Y3 * Y3 + Z3 * Z3) 

IAM = arccos(Z3/MAG) 

11. Calculate the vis-viva parameter. 

QQM = R1 * VELM * VELM/|i m 

12. Calculate eccentric anomaly. 

COSAEX = (EEM + cos(0))/(l + EEM * cos(©)) 

AEX = arccos(COSAEX) 

13. Calculate semi-major axis. 

AAX = Rl/(2 - QQM) 

14. Store old "flyby" time of flight. 

TIEMMS = TIEMM 

15. Calculate the time of perigee passage. 

A. If the orbit is elliptical: 

TIMX = (AAX A 3/p tt ) * (AEX - EEM * sin(AEX)). 

B. If the orbit is hyperbolic: 

TIMX = (-AAX A 3/pJ * (EEM * sinh(EEM) 
sinh(EEM))). 

16. Calculate time of flight for "flyby". 

TIEMM = TIMX * 2 

17. Calculate longitude of nodes. 

ANODE = OMEGM * TIEMM 
IF ((ALAT*MOOD) >= 0) then ANODE = ANODE + PI 

18. Calculate longitude angle between nodes and the SOI projection. 

AVM = ALON - ANODE 


log(cosh(EEM) + 



19. Calculate angle between nodes and SOI projection. 

COSPHIM = cos(ALAT) * (cos(ALODE) * cos(ALON) + sin(ANODE) * 
sin(ALON)) 

APHIM = arccos(COSPHIM) 

20. Calculate azimuth angle for R2 (at SOI). 

SINAZM = cos(AIM)/cos(ALAT) 

COSAZM = sin(AIM)*cos(APHJM)/cos(ALAT) 

AZMM = arctan2(SINAZM, COSAZM) 

21. Calculate position for LI (X-Y axes are fixed at Moon) 

ALATL1 = 0 

ALONL1 = arctan2( Y1, XI) 

22. Calculate angle between nodes and LI. 

COSPHIL1 = cos(ALATLl) * (cos(ALODE) * cos(ALONLl) + sin(ANODE) * 
sin(ALOLlN)) 

APHIL1 = arccos(COSPHILl) 

23. Calculate azimuth angle for R1 (at LI). 

SINAZM = cos(AIM)/cos(ALATLl) 

COSAZM = sin(AIM)*cos(APHILl)/cos(ALATLl) 

AZMLP1 = arctan2(SINAZM, COSAZM) 

24. If new "flyby" flight time is not very close to the saved value go back to step 2. 

25. Set latitude/longitude position of LI (X-Y axes are rotating) 

ALATLP1 = 0° 

ALONLP1 = 180° 

26. Call Veltrans to establish the velocity components of LI in Lunar-based rectangular 
components. 

27. Set Y-component of velocity for LI (X- and Z- components are 0). 

VYCOR = -OMEGM * Rl/3600. 

28. Calculate the AV needed to correct the velocity vector at LI to that of Li’s orbit 
(outbound) or the "flyby" to the SOI (return). 

DELYLP1 = (VXLP 2 + (VYLP1 - VYCOR) 2 + VZLP1 2 ) 

29. Sign of flight path angle must be changed for SOI calculations. 

GAMAM = -GAMAM 




Appendix G. Subroutine GAMACALC 


The subroutine GAMACALC receives the parameters orbital perigee radius (RPX), orbital 
velocity at LP (VV), orbital radius at LP (RRX), and the gravitational parameter of the body 
being orbited (CMUX). It calculates and returns time of flight (TEMX), the cosine of the flight 
path angle (COSGAMX), velocity at periapses (VPX), circular velocity at periapses (VCIRX), 
true anomaly (THETAX), and an indicator describing whether the orbit is elliptical or hyperbolic 
(TRAJ$). 

1. Initialize indicators to presume an elliptical orbit. (IHYPER, TRAJ$). 

2. Calculate the vis-viva parameter 

QQX = (RRX * VELX A 2) / CMUX. 

3. If the orbit is just barely hyperbolic (QQX is within one-millionth of 2), force the 
calculation to consider it elliptical (reduce QQX to just under 2). 

4. If the orbit is still hyperbolic, reset the indicators to show this. Print a message on the 
screen announcing a hyperbolic orbit. 

5. Calculate the semi-major axis of the orbit 

AAX = RRX / (2-QQX). 

6. If the semi-major axis is very large (greater than 10 A 12) or very small (less than -10 A 12), 
the orbit is trapped near a parabolic trajectory. Make it hyperbolic: 

AAX - -10 A 12. 

7. Calculate the orbit eccentricity 

EEX = 1 - (RPX/AAX). 

8. Calculate the semi-latus rectum 

PPX = AAX * (1 - EEX). 

9. Calculate the flight path angle 

a. The angular momentum of the orbit at perigee is 

HP = RPX * VELPERIGEE * cos(GAMAX). 

But at perigee, GAMAX is zero, so 
HP = RPX * VELPERIGEE. 

b. At LP, angular momentum is 

HX = RRX * VELX * cos(GAMAX). 

c. Angular momentum is constant along a given orbit, so 

HP = HX 

RPX * VELPERIGEE = RRX * VELX * cos(GAMAX) 

GAMAX = arccos((RPX * VELPERIGEE) / (RRX * VELX)) 




d. The velocity at perigee is 

VELPERIGEE = sqr((CMUX * RAPOGEE) / (AAX * RPX)). 

Since RAPOGEE / AAX = (1 + EEX), then 

VELPERIGEE = sqr (CMUX * (1 + EEX) / RPX) or 
RPX * VELPERIGEE = sqr(RPX * (1 + EEX) * CMUX). 

e. Substituting (d) into (c) above yields 

GAMAX = arccos((sqr(RPX * (1 + EEX) * CMUX) / (RRX * VELX)) 

f. Substituting from (2) above yields 

GAMAX = arccos(sqr((RPX * (1 + EEX)) / (RRX * QQX))). 

10. Calculate the perigee velocity 

VPX = sqr(CMUX * (1 - EEX) / RPX). 

11. Calculate the circular velocity at perigee 

VCIRX = sqr(CMUX / RPX). 

12. Calculate the true anomaly 

THETAX = arccos(((PPX / RRX) -1) / EEX) 

13. Calculate the eccentric anomaly 

AEX = arccos((EEX + cos(THETAX)) / (1 + EEX * cos(THETAX))). 

14. Compare the orbital radius at LP (RRX) to the apogee radius (AAX * (1 + EEX)). If the 
orbital radius at LP is greater than the apogee radius, print a message on the screen 
indicating the difference in nautical miles. Also display the following: 

a. LP latitude (ALAT) 

b. LP longitude (ALON) 

c. vis-viva parameter (QQX) 

d. semi-major axis (AAX) 

e. eccentricity (EEX) 

f. flight path angle (GAMAX) 

g. perigee velocity (VPX). 

15. Calculate the time of flight. 

a. If the orbit is elliptical: 

TTMX = sqr(AAX A 3/CMUX) * (AEX - EEX * sin(AEX)). 

b. If the orbit is hyperbolic: 

TIMX = sqr(-AAX A 3/CMUX)*(EEX*sinh(EEX) - log(cosh(EEX) + 
sinh(EEX))). 




Appendix H. Subroutine POSVELMO 


The subroutine POSVELMO receives the parameter TIM (number of Julian centuries from the 
year 2000) and returns the moon’s position and velocity at that time. Specifically, it returns the 
moon’s distance from the Earth’s center, in Earth radii (RRM); velocity in the x-direction (along 
the Earth-Moon line), in feet per second (XDLO); and velocity in the y-direction (direction of 
the moon’s orbit), in feet per second (YDLO). 

1. Calculate a fraction of time that represents half a second. 

DELT = 0.5 seconds / (36525 days/century * 24 hrs/day * 3600 seconds/hr) 

= 1.58440E-10 centuries. 

2. Call the subroutine MOON with the parameter (TIM minus DELTA) to determine the 
moon’s distance in Earth radii at half a second before TIM. This distance is RM1. 

3. Call the subroutine MOON with the parameter (TIM plus DELTA) to determine the 
moon’s distance in Earth radii at half a second after TIM. This distance is RM2. 

4. Calculate the velocity of the moon in the -x (radial) direction. 

XDLO = [(RM2 Earth radii - RM1 Earth radii) / (1 sec)] * 20,925,741 ft/radii. 

5. Calculate the average radius of the Lunar orbit during the one second centered on TIM. 

RRM = (RM1 + RM2) / 2. 

6. Determine the radius of the Lunar orbit from the Earth-Moon barycenter. 

RRMB = RRM - 0.7412789 Earth radii. 

7. Calculate the moon’s velocity in the direction of its orbit. 

a. Moon’s apogee (Apo) = 62.83308 Earth radii. 

b. Moon’s perigee (Per) = 55.68264 Earth radii. 

c. Eccentricity (e) = (Apo - Per) / (Apo + Per) = 0.06033. 

d. Semi-latus rectum (p) = Apo(l -e A 2) = 62.60439 Earth radii 

= 1,310,038,967 feet. 

e. Earth’s gravitational parameter (mu) = 1.407646822E+16 ft A 3/sec A 2. 

f. Angular momentum (h) = sqr(mu * p) 

= 4.29427E+12 ft A 2/sec * (1 earth radii / 20,925,672.57 ft) 

= 205,215.4 ft*Earth radii/second. 


Y-velocity (YDLO) = h / RRMB = 205,215.4/RRMB. 




Appendix I. Subroutine MOON 


The subroutine MOON receives the parameter T (number of Julian centuries from the year 2000) 
and returns the approximate location of the moon in geocentric coordinates at that time. 
Specifically, it returns the right ascension of the moon (RAM), declination of the moon 
(DECM), and distance to the moon in Earth radii (RM). The formulae are from The Astronomi¬ 
cal Almanac of the Year 1984 . page D46. All degrees are converted to radians with the 
conversion factor C = tc/180. 

1. Calculate the ecliptic coordinates of the moon. 

a. Moon’s ecliptic longitude 

LAM = 218°.32 + 481267°.833T 

+ 6°.29 * sin(134°.9 + 477198°.85T) 

-1°.27 * sin(259°.2 - 413335°.38T) 

+ 0°.66 * sin(235°.7 + 890534°.23T) 

+ 0°.21 * sin(269°.9 + 954397°.70T) 

- 0°.19 * sin(357°.5 + 35999°.05T) 

- 0°.ll * sin(186°.6 + 966404°.05T) 

b. Moon’s ecliptic latitude 

BETA = 5°.13 * sin( 93°.3 + 483202°.03T) 

+ 0°.28 * sin(228°.2 + 960400°.87T) 

- 0°.28 * sin(318°.3 + 6003°. 18T) 

- 0°.17 * sin(217°.6 - 407332°.20T) 

c. Horizontal parallax 

PEE = 0°.9508 

+ 0°.0518 * cos(134°.9 + 477198°.85T) 

+ 0°.0095 * cos(259°.2 - 413335°.38T) 

+ 0°.0078 * cos(235°.7 + 890534°.23T) 

+ 0°.0028 * cos(269°.9 + 954397°.70T) 

d. Semi-diameter of moon’s orbit 

SD = 0.2725 * PIE 

e. Distance to the moon in Earth radii 

RM = 1 / sin(PIE) 

2. Form the geocentric direction cosines to rotate into geocentric coordinates. 

a. 1 = cos(BETA)cos(LAM) 

b. m = 0.9175*cos(BETA)sin(LAM) - 0.3978*sin(BETA) 

c. n = 0.3978*cos(BETA)sin(LAM) + 0.9175*sin(BETA) 

where 1 = cos(DECM)cos(RAM), m = cos(DECM)sin(RAM), n = SIN(DECM). 





3. 


Then: 


a. RAM = arctan(m/l) [right ascension] 

b. DECM = arcsin(n) [declination] 

The errors will rarely exceed 0.2 Earth radii in distance (RM), 0.3° in right ascension (RAM), 
and 0.2° in declination. 




Appendix J. Subroutine VELTRANS 


The subroutine VELTRANS converts an orbital velocity vector into rectangular coordinates (see 
figure K-5). Parameters received by the subroutine are velocity (VEL), flight path angle 
(GAMA), azimuth (AZM), latitude above the Earth-Moon plane (ALAT), and longitude from 
the Earth-Moon line (ALON). A set of intermediate calculations is performed to express the 
velocity vector in terms of a radial component, a latitudinal component, and a longitudinal 
component. Each of these three components is further resolved into x-, y-, and z-components. 
Finally, all three x-components, all three y-components, and all three z-components are summed 
to provide the total x-, y-, and z-components of velocity. 

1. Conversion of velocity vector into spherical coordinates. 

From the geometry, the radial component of velocity, R, is calculated to be 
VEL * sin(GAMA). (See figure K-6). 

The component along the orbital path, R , is 
VEL * cos(GAMA). 

This orbital path component of velocity is further resolved into a latitude component, 

LAT, and a longitude component, LON (see figure K-7). Again, from the geometry, 

LON = R * sin(AZM) and 
LAT = R * cos(AZM). 

2. Conversion of spherical coordinates into rectangular coordinates. 

a. Conversion of radial component into rectangular coordinates. 

Refer to figure K-8. The projection of R onto the x-y plane is 
R * cos(LAT). 

The negative x-component of this is 
R * cos(LAT) * cos(LON) 
so the x-component, X, is 

-R * cos(LAT) * cos(LON). 

The negative y-component of R is 
R * cos(LAT) * sin(LON) 
so the y-component, Y, is 

-R * cos(LAT) * sin(LON). 

The z-component, Z, is 
R * sin(LAT). 

b. Conversion of latitude component into rectangular coordinates. 

Refer to figures K-9 and K-10. LAT is perpendicular to the radial vector, R. A 
line in the z-direction that meets the tip of LAT and intersects the ^pjdius vector 
produces the angles a and b, where 
b = 90 - LAT and 
a + b + 90 = 180. Therefore, 
a = LAT. 

From the geometry, the z-component of LAT, ZLAT, is 
LAT * cos(LAT). 




The projection of LAT onto the x-y plane is 
LAT * sin(LAT). (see figures K-11). 

The x-component of this, XLAT, is 
LAT * sin(LAT) * cos(LON). 

The y-component of this, YLAT, is 
LAT * sin(LAT) * sin(LON). 

c. Conversion of longitude component into rectangular coordinates. 

Refer to figures K-12 and K-13. LON is always parallel to the x-y plane, so the 
z-component of LON, ZLON, is always zero. Using the same proof as in (b) 
above, it can be seen that the angle between LON and the y-component of LON is 
equal to LON. From the geometry, the x-component of LON, XLON, is 
LON * sin(LON). 

The negative y-component of LON is 
LON * cos(LON), 
so the y-component, YLON, is 
-LON * cos(LON). 

Sum of the rectangular coordinates. 

All of the x-, y-, and z-components are summed to provide the complete rectangular 

coordinates of die velocity vector. 

VXX = X +XLAT + XLON 
VYX = Y + YLAT + YLON 
VZX = Z + ZLAT + ZLON. 




The projection of LAT onto the x-y plane is 
LAT * sin(LAT). (see figures K-l 1). 

The x-component of this, XLAT, is 
LAT * sin(LAT) * cos(LON). 

The y-component of this, YLAT, is 
LAT * sin(LAT) * sin(LON). 

c. Conversion of longitude component into rectangular coordinates. 

Refer to figures K-l2 and K-l3. LON is always parallel to the x-y plane, so the 
z-component of LON, ZLON, is always zero. Using the same proof as in (b) 
above, it can be seen that the angle between LON and the y-component of LON is 
equal to LON. From the geometry, the x-component of LON, XLON, is 
LON * sin(LON). 

The negative y-component of LON is 
LON * cos(LON), 
so the y-component, YLON, is 
-LON * cos(LON). 

Sum of the rectangular coordinates. 

All of the x-, y-, and z-components are summed to provide the complete rectangular 

coordinates of the velocity vector. 

VXX = X + XLAT + XLON 
VYX = Y + YLAT + YLON 
VZX = Z + ZLAT + ZLON. 




Appendix K. Figures 







Figure K-3 


















Figure K-12 





